General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



NASA Technical Memorandum 82876 


(NASA-Ttt-d2B7t>) xftPLICIT aARCUING SOLUTION 
Of COttPa£SSIt>LL VISCOUS SUBSONIC FLOW IN 
PLANAR AND A aIS V ftHETRIC DUCTS Ph.D. Tbesis 
.NASA) z22 p riC AlO/flF AU1 CSCL 2JD 

Gi/ia 


NPJ- 103B9 


Unclas 

3839d 


Implicit Marching Solution of 
Compressible Viscous Subsonic 
Flow in Planar and 
Axisymmetric Ducts 


y 

Charles E. Towne 
Lewis Research Center 
Cleveland , Ohio 

and 



Joe D. Hoffman 
Purdue University 
West Lafayette, Indiana 


September 1982 





1 


TABLE Or CONTENTS 

Page 

PRINCIPAL NOTATION 

SUMMARY » . lx 

SECTION 

1. INTROUUCTION 1 

1.1 GENERAL BACKGROUND 1 

1.2 EARLIER WORK 2 

1.3 PRESENT ANALYSIS 9 

2. ANALYSIS 11 

2.1 COOROINAlt SYSTEM 11 

2.2 GOVERNING EQUATIONS 13 

2.2.1 Equations of Motion 13 

2.2.2 Boundary Conditions 21 

2.2.3 Initial Profiles 2s 

2.2.4 Molecular Transport Properties ..... 28 

2.3 STREAMWlSE PRESSURE GRADIEN1 29 

2.3.1 Viscous Pressure Correction 2» 

2.3.2 Imposed Pressure field 31 

2.4 TURBULENCE MODEL 33 

2.4.1 Outer Region Model 34 

2.4.2 Cebeci-Smith Inner Region Model 3s 

2.4.3 McDonald-Camarata inner Region Mooel 37 

3. NUMERICAL METHOD 41 

3.1 BASIC DIFFERENCING PROCEDURE 4l 

3.1.1 Computational Mesh 41 

3.1.2 Difference Formulas 43 

3.2 L IhEARiZAl ION PROCEDURE 44 

3.3 BOUNDARY CONDITIONS bl 

3.3.1 Derivatives at Boundaries bl 

3.3.2 Symmetry Line b2 

3.3.3 Solia Surfaces b3 

3.3.4 Wall Functions s4 

3.4 TRI-DIAGONAL MATRIX INVERSION ALGORITHM bb 

3.5 VISCOUS PRESSURE CORRECTION SB 

3.5.1 Streamwise Momentum Equation bd 

3.5.2 Total Mass Flow Rate Equation S9 

3. b. 3 Solution Algorithm 61 


11 


4. RESULTS 6S 

4.1 LAMINAR DEVELOPING PIPE FLOW 65 

4.1.1 Grid and Initial Conditions 66 

4.1.2 Results 68 

4.2 JEFFERY-HAMEL FLOW 81 

4.2.1 Grid and Initial Conditions 81 

4.2.2 Results 83 

4.3 TURBULENT DEVELOPING PIPE FLOW 88 

4.3.1 Grid and Initial Conditions. . 88 

4.3.2 Results 90 

4.4 S-SHAPED DUCT FLOW 100 

4.4.1 Geometry and Computational Coordinates 100 

4.4.2 Imposed Pressure Field 102 

4.4.3 Initial Conditions 103 

4.4.4 Results 104 

4.5 DIFFUSER FLOW 114 

4.5.1 Geometry and Computational Coordinates 114 

4.5.2 Imposed Pressure Field 116 

4.5.3 Initial Conditions 117 

4.5.4 Results 118 

5. CONCLUSIONS 132 

LIST OF REFERENCES 134 

APPENDICES 

A. DERIVATION OF GOVERNING EQUATIONS 140 

B. ORDER-OF -MAGNITUDE ANALYSIS 152 

C. POTENTIAL FLOW COMPRESSIBILITY CORRECTION 156 

D. COUPLED DIFFERENCE EQUATIONS 158 

E. DENSITY BOUNDARY CONDITION AT A WALL 180 

F. COEFFICIENT MATRICES AND SOURCE TERMS 183 

G. DIFFERENCE EQUATIONS FOR VISCOUS PRESSURE CORRECTION 198 

N. JEFFERY-HAMEL FLOW STREAMWISE PRESSURE GRADIENT 209 



PRINCIPAL NOTATION 


Coefficient of velocity in uncoupled streamwise 
momentum equation 

Area 

Coefficient matrix in coupled equations 

Coefficient of velocity in uncoupled streamwise 
momentum equation 

Coefficient matrix in coupled equations 

Coefficient of velocity in uncoupled streamwise 
momentum equation 

Skin friction coefficient 

Specific heat at constant pressure 

Pressure coefficient 

Specific heat at constant volume 

Coefficient in coupled continuity equation 

Coefficient in coupled energy equation 

Coefficient in cross-flow momentum equation 
evaluated at a wall 

Coefficient in coupled streamwise momentum 
equation 

Coefficient in coupled cross-flow momentum 
equation 

Coefficient matrix in coupled equations 

Coefficient of viscous pressure correction in 
uncoupled streamwise momentum equation 

Pipe diameter 

Parameter controlling grid packing 


1v 


e 


e mn 

f 

h 

hi,h 2 ,h 3 

i 

j 

J 

k 

t 

*-e 

L r 

• 

m 

M 

Mi 

C 

P 

Pu 

Pret 


P' 

P 

Pr 

r 

r F 


Coefficient of velocity in total mass flow rate 
equation 

Strain rate tensor (m and n - 1* 2, 3) 

Coefficient of viscous pressure correction in 
total mass flow r*te equation 

Static enthalpy 

Metric scale coefficients 

Grid index in streamwise U) direction 

Grid index in cross-stream (n) direction 

Value of j o', upper wall 

Thermal conductivity coefficient 

Mixing K igth 

entrance length in developing pipe flow 
Reference length 
Mass flow rate 
Mach number 

Nominal inlet Mach number 
Symbol for "order of" 

Static pressure 
Total pressure 

Reference pressure in definition of pressure 
coefficient 

Viscous static pressure correction 
Imposed static pressure 
Prandtl number 

Heat flux vector (in ■ 1, 2, 3) 

Absolute Cartesian or polar coordinate 
Recovery factor 


V 


R 

Reg 


Re r 

Re R 


S M 

$x 


Sw 

s x 




s 

t 

T 

t aw 

V 

T + 

u 


u 0 


u 


T 


+ 


U 


* 


u 


Gas constant; Pipe radius 

Reynolds number based on pipe diameter and 
average velocity 

Reference Reynolds number 

Reynolds number based on pipe radius and 
average velocity 

Source term in total mass flow rate equation 

Source term in uncoupled streamwise momentum 
equation 

Source term in coupled continuity equation 

Source term in coupled energy equation 

Source term in cross-flow momentum equation 
evaluated at a wall 

Source term in coupled streamwise momentum 
equation 

Source term in coupled cross-flow momentum 
equation 

Source term vector in coupled equations 
Time 

Static temperature 
Total temperature 
Adiabatic wall temperature 
Reference temperature 
Law-of-the-wal 1 temperature 
Streamwise velocity 
Average streamwise velocity in pipe 
friction velocity 

Law-of-the-wal 1 streamwise velocity 

Streamwise velocity from preliminary ...arcning 
step 


Vi 


U r 

v 

Vp 

w 

X 

y 

y* 

Y 


z 



6 


i 

Y 


& 


* 


e i 

c o 

n 

e 


K 

X 


V 

v 

C 

p 

p r 


Reference velocity 

Cross-flow, or normal, velocity 

Potential flow velocity 

Implicit weighting factor 

Solution vector in coupled equations 

Distance from wall 

Law-of-the-wall coordinate 

Transformed (equally spaced) cross-stream 
coordinate 

Absolute Cartesian or polar coordinate 

Ratio of specific heats, c p /c v 

boundary layer thickness 

Kronecker delta (m and n - 1, 2, 3) 

Difference operator for first derivative in 
Y-direction 

Difference operator for second derivative in 
Y-direction 

Boundary layer displacement thickness 
Inner region eddy viscosity 
Outer region eddy viscosity 
Computational cross-stream coordinate 
boundary layer momentum thickness 
Von Karman constant 

Ratio of successive streamwise step sizes 
Viscosity coefficient 
Kinematic viscosity coefficient 
Computational streamwise cooroinate 
Static density 
Reference density 


vl 1 


T mn 

v 

Subscripts 

AVE 

C-8QDY 

COWL 

e 

E 

F-D 

i 

j 

L 

r 

T 

W 

n 

E 

S uperscript 

T 


Stress tensor (m and n » i, 3) 

Viscous dissipation 

Average at inlet station 
Along centerbody 
Along cowl 

At edge of boundary layer 
Effective (laminar plus turbulent) 
Fully-developed 

Grid location in streamwise (0 direction 

Grid location in cross-stream (n) direction 

Laminar 

Reference 

Turbulent 

Wall 

Partial differentiation in n-direction 
Partial differentiation in E-direction 

Vector transpose 


ix 


IMPLICIT MARCHING SOLUTION OF COMPRESSIBLE VISCOUS 
SUBSONIC FLOW IN PLANAR AND AXISYMMETRIC DUCTS 

Charles E. Towne 

National Aeronautics and Space Administration 
Lewis Research Center 

Cleveland, Ohio 

and 

Joe D. Hoffman 
Purdue University 
West Lafayette, Indiana 

SUMMARY 

A new streamwise marching procedure has been developed and coded 
for compressible viscous subsonic flow in planar or axisymmetric ducts 
with or without centerbodies. The continuity, streamwise momentum, 
cross-flow momentum, and energy equations are written in generalized 
orthogonal curvilinear coordinates. To allow the use of a marching 
procedure, second derivatives in the streamwise direction are neglec- 
ted, and the pressure in the streamwise momentum equation is written 
as the sum of a known two-dimensional imposed pressure field and an 
unknown one-dimensional viscous correction. For turbulent flow, the 
Reynolds stress and turbulent heat flux terms are modeled using two 
different two-layer eddy viscosity turbulence models. 

Prior co each main marching step, a preliminary marching step is 
taken irt which the integral mass flow rate equation and an uncoupled 
form of the streamwise momentum equation are solved simultaneously to 
obtain the viscous pressure correction. During the main marching step 
the four governing equations are solved simultaneously as a coupled 
system using an implicit finite-difference method, with the viscous 
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pressure correction treated as a source term. The equations are 
linearized to second-order accuracy without using iteration by expand- 
ing each unknown term in a Taylor series in the streamwise direction. 

Results are presented for developing laminar flow in a circular 
pipe, laminar flow in a two-dimensional converging channel (Jeffery- 
Hamel flow), developing turbulent flow in a circular pipe, turbulent 
flow in a two-dimensional S-duct, and turbulent flow in a typical sub- 
sonic diffuser for a supersonic inl^t. For all test cases, the results 
are compared with data and/or exact solutions. 

The computed results agree very well with the data and exact so- 
lutions for th laminar cases studied. The results for the turbulent 
cases also agree well with the data, although not quite as well as for 
laminar flow. The viscous pressure correction properly accounts for 
the effect of viscous blockage on the streamwise pressure gradient. 
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SECTION 1 
INTRODUCTION 


1.1 GENERAL BACKGROUND 


The computation of viscous subsonic internal flow has many engi- 
neering applications. The design or analysis of a subsonic jet engine 
system, for example, requires the computation of flow through the inlet 
and nozzle. Even a supersonic inlet requires a subsonic solution in 
the diffuser. 

The traditional method for computing subsonic flow in a duct is to 
first use an inviscid analysis to get a pressure field in the main 
core flow region, and then to use a boundary layer analysis to account 
fo , “ the viscous effects near the walls. If the boundary layers grow 
large enough to have a blockage effect on the core flow, however, some 
type of patching or iteration between the inviscid and boundary layer 
analyses must be done. If the boundary layers grow further, so that 
they merge ar.d fill the duct with viscous flow, this method breaks 
down completely. 

One approach that can be used in this situation is to solve the 
complete Navier-Stokes equations. Since these equations are elliptic, 
they are usually solved by time-dependent or relaxation techniques, 
requiring large amounts of computer time and storage. 
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For many flows, however, a complete Navier- Stokes solution is not 
necessary. If the flow conditions and geometry are such that elliptic 
effects in one of the coordinate directions can either be ignored or 
assumed known a priori, a marching solution procedure can be used. 

This requires that viscous and thermal diffusion in the streamwise 
direction be neglected. In addition, special treatment of the stream- 
wise pressure gradient is normally required, since the pressure field 
in subsonic flow is inherently elliptic. Physically, these assumptions 
imply a relatively high Reynold's number flow without abrupt changes 
in geometry and without flow separation. 

1.2 EARLIER WORK 


All previously published marching methods for viscous internal 
flow, of course, neglect streamwise diffusion. However, they differ 
considerably in several details, including the form of the equations 
being solved, the degree of approximations and/or simplifications in- 
volved, the treatment of the streamwise pressure gradient, the degree 
of coupling between equations, the treatment of compressibility, and 
the numerical method used. 

An early effort in this area was that of Anderson (Refs. 1-3). He 
developed a method for analyzing compre'sible axisymmetric flow with 
or without swirl by assuming that the cross-flow velocity v is of 
order 6 and then neglecting all C (6 ) terms. In addition, sev- 
eral other terms are neglected that involve either derivatives of 
metric scale coefficients or streamwise derivatives of v. An inviscid 
cross-flow momentum equation is used, with all of the convective terms 
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neglected (for zero swirl) except the one involving the derivative of 
the metric scale coefficient in the cross-flow direction. Thus, the 
cross-flow pressure grad int is balanced by the curvature of the coor- 
dinate lines parallel to the wall. The coordinate l'nes are given by 
the streamlines and potential lines from a planar potential flow solu- 
tion. In his orininal work (Ref. 1), Anderson uncoupled the equations 
and solved them using an explicit numerical method. He later modified 
the analysis by expanding the dependent variables in a Taylor series 
in the marching direction and solving the resulting equations impli- 
citly (Refs. 2-3). In this formulation, the equation of state and the 
definitions of the two shear stresses, heat flux, and entropy (the 
energy equation dependent variable) are treated as independent equa- 
tions. The implicit solution is thus accomplished by the inversion of 
a block tri-diagonal matrix with 10x10 sub-matrices. This method has 
been widely applied to flows in subsonic diffusers (Refs. 4-7). 

In Anderson's method, because of the simplifications made to the 
governing equations, the st.reamwise pressure gradient p^ can be 
treated like any ether term. However, as demonstrated by several other 
authors (Refs. 0-11), when a more complete form of the cross-flow mo- 
mentum equation is solved, special treatment of the streamwise pressure 
gradient is required. 

In the three-dimensional method of Patankar and Spalding (Ref. 9), 
the streamwise pressure gradient p^ is assumed to be a function 
of streamwise distance only. Thus, the pressure gradient in the 
streamwise momentum equation is uncoupled from those in the cross-flow 
momentum equations. To advance the solution one marching station using 
the Patankar-Spalding method, the procedure is as follows: 
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1. Guess a value for p^U). 1° practice the value for the 

preceding step would normally be used. 

2. Solve the streamwise momentum equation, uncoupled from the rest 
of the equations, for the streamwise velocity u using an implicit 
numerical method. Since the continuity equation has not yet been used, 
the computed mass flow rate will, in general, be incorrect. 

3. Correct p^., using an integral equation derived from an 
approximate inviscid form of the streamwise momentum equation, to get 
the correct total mass flow rate. 

4. Correct u to be consistent with step 3. 

b. Guess a pressure distribution in the new cross-flow plane. 

Again, in practice the values at the preceding station would be used. 

6. Solve the cross-flow momentum equations, uncoupled, for the 
cross-flow velocities v and w. 

7. Correct the cross-flow pressures by solving a two-dimensional 
Poisson equation derived from the continuity equation and approximate 
forms of the inv;scid cros^-flow momentum equations. 

8. Correct v and w to be consistent with step 7. 

9. Solve the energy equation, uncoupled, for the temperature dis- 
tribution at the new station. 

If the flow is compressible, the density would presumably be updated 
after step 9, although this is not clear in Reference 9. Multiple-sweep 
variations of this method, terued "partially-parabolic" procedures, 
have also been used (Refs. 12-14). However, by sweeping the entire 
flow fieid several times, much of the advantage of a marching solution 
over a complete Navier- Stokes solution is lost. The Patankar-Spalding 
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method has been applied to straight pipes with circular, square, and 
polar cross-sections (Refs. 9, 11, 15-16), to curved circular pipes 
(Refs. 17-19), and, with iteration at each station to reduce lineari- 
zation errors, to a 0-shaped transitioning duct (Ref. 15). 

The method of Briley (Ref. 20) is similar in concept to that of 
Patankar and Spalding, but is somewhat more rigorous and differs con- 
siderably in detail. He splits the pressure into an inviscid contri- 
bution P and a viscous correction p*. The inviscid contribution is 
assumed to be known a priori, and is used to bring elliptic effects of 
the geometry into the solution. The pressure correction term in the 
streanwise momentum equation is uncoupled from the correction terms in 
the cross-flow momentum equations in the manner suggested by Patankar 
and Spalding. The method for computing this pressure, however, is 
different. Briley's procedure for advancing the solution one step is 
as follows: 

1. Guess a value for p^U). Again, in practice it would be 
lagged one step. 

2. Solve the streamwise momentum equation, uncoupled, by an impli- 
cit method. 

3. Correct p^, using secant iteration on steps 1 and 2, to get the 
correct total mass flow rate. 

4. Solve the cross-flow momentum equations (uncoupled), using a 
lagged pressure field, to get predicted cross-flow velocities v p 
and Wp. 

5. Solve a two-dimensional Poisson equation for where ? 

is a velocity potential for the cross-flow velocity corrections v^ 
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and This Poisson equation is derived from the continuity equa- 

tion by assuming the cross-flow velocity corrections are irrotational. 

b. Compute Vq and w^ from *. Note that with this formulation 
v c and w^ will have non-zero values at the wall, and therefore the 
final cross-flow velocities will not satisfy the no-slip boundary con- 
dition. 

7. Correct the pressure field used in the cross-flow momentum equa- 
tions by solving a two-dimensional Poisson equation derived from the 
cross-flow momentum equations. 

8. Solve the energy equation for the temperature at the new sta- 
tion. 

As presented in Reference 20, Briley's method is for incompressible 
flow, so the requirement for updating the density does not arise. 

This method was later improved by Briley and McDonald (Ref. 21) 
and by Levy, McDonald, Briley, and Kreskovsky (Ref. 22). In the newer 
method, the velocity is split into a primary (streamwise) component 
and a secondary (cross-flow) component. The secondary flow component 
is split further into irrotational and rotational contributions, which 
are presumed derivable from a secondary flow velocity potential and 
stream function, respectively. The pressure is split as in Briley's 
original method. To march one step with this method, the procedure is 
as follows: 

1. Guess p^U) . 

2. Solve the streamwise momentum equation, the energy equation, 
and an approximate equation of state for the primary velocity u, the 
static enthalpy h, and the density p. In many applications, the 
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total enthalpy is assumed constant, which eliminates the requirement 
for solving the energy equation. 

3. Correct p^, using secant iteration on steps i and 2, to get the 
correct total mass flow rate. 

4. Solve a two-dimensional Poisson equation for *, where t 

is the velocity potential for the irrotational component of the cross- 
flow velocity. This Poisson equation is derived from the continuity 
equation. 

5. Solve an approximate transport equation for q and a two- 
dimensional Poisson equation for * as a coupled system, where 

a is the streamwise vorticity and * is the stream function for 
the rotational component of the cross-flow velocity. The streamwise 
vorticity equation can be derived from the cross-flow momentum equa- 
tions. The Poisson equation for * is simply the definition of 
streamwise vorticity written in terms of the secondary flow stream 
function. 

b. Compute the cross-flow velocities from * and *. 

With this formulation, no explicit correction is made to the pressure 
in the cross-flow momentum equations, since it does not appear in the 
vorticity transport equation. At the end of a calculation one could 
solve a three-dimensional Poisson equation for the pressure field since 
the velocities are known. In practice, however, this has not been 
done. The method has been applied to straight circular pipes 
(Ref. 22), curved ducts of rectangular cross-section (Refs. 21,23), a 
turbine blade passage (Ref. 21), and mixer nozzles (Refs. 24-26). 
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A method has been developed by Dodge (Ref. 27) for three-dimen- 
siunal incompressible flow that is quite different from the preceding 
methods. He splits the velocity into inviscid and viscous components 
V. and V v , respectively, and assumes that the inviscid component 
is derivable from a velocity potential f. That is, 

» * » > 

V « V- ♦ V « + V 

i v T v 

This expression is substituted into the three momentum equations. It 
is then assumed that the pressure gradients in all three directions 
are balanced by the purely potential terms. Thus, the pressure 

gradients are eliminated from the momentum equations, which can then 

> 

be solved for V y by a marching technique if 9 is known. When the 
assumed velocity split is substituted into the continuity equation, the 
result is a three-dimensional Poisson equation for 9 if is 
known. The numerical procedure is thus: 

1. Guess a three-dimensional potential fiela (or, equivalently, a 
pressure field) . 

2. Solve the three momentum equations throughout the flow field for 

* 

V v using a marching procedure. It is not clear in Reference 27 
whether an explicit or implicit method is used. 

3. Solve the three-dimensional Poisson equation, using the results 
for V v from step 2, to get a new potential field. 

4. Iterate steps 2 and 3 until convergence. 

Note that this method involves multiple sweeps of the computational 
domain. As with the "partial ly-parabol ic" variations of the Patankar- 
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Spalding method, if more than a few sweeps are needed, a complete 
Navier-Stokes solution might be preferred. 


1.3 PRESENT WORK 


In the present study, a new streamwise marching procedure was de- 
veloped and coded for compressible viscous subsonic flow in planar or 
axisymmetric ducts with or without centerbodies. As in the previously 
described methods, second derivatives in the marching direction are 
neglected. The remaining second-order viscous terms are retained, 
however, since they can be included within the framework of a marching 
procedure. 

Following Briley (Ref. 20), the pressure in the streamwise momen- 
tum equation is written as the sum of a known two-dimensional inviscid 
contribution P and an unknown one-dimensional viscous correction p*. 
However, this form for the pressure is used only in the streamwise 
momentum equation, not in the cross-flow momentum equation. This is 
sufficient to remove the streamwise elliptic character of the pressure 
field and allow the use of a marching solution procedure. The pressure 
in the cross-flow momentum equation is simply treated as an unknown 
dependent variable, thus eliminating the requirement for a separate 
correction procedure for the cross-flow pressure field. Since elliptic 
effects due to geometry are included in the known inviscid pressure 
field, only one sweep through the flow field is required. 
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In order to rigorously model the relationships between dependent 
variables* the four governing equations (continuity, streamwise momen- 
tum, cross-flow momentum, and energy) are solved simultaneously at 
each station as a coupled system. Since the continuity equation is 
one of the equations in the system, however, conservation of total 
mass flow rate cannot be used during or after a marching step to com- 
pute the pressure correction p*. Therefore, a preliminary marching 
step is taken using the integral mass flow rate equation and an un- 
coupled form of the streamwise momentum equation. Instead of iterating 
to find the value of p' that gives the correct mass flow rate, these 
two equations are solved simultaneously as a coupled system to get p' 
directly. Then, with p' known, the four governing equations are 
solved simultaneously with p^ treated as a source term. 

The equations are solved by an implicit finite-difference proce- 
dure. A weighting factor is used to allow the degree of implicitness 
to vary from Crank-Nicholson type to fully implicit. The equations are 
linearized by expanding each unknown term in a Taylor series in the 
streamwise direction. This allows the inherent nonlinearities of the 
governing differential equations to be modeled with second-order accu- 
racy without using iteration. Cross-flow derivatives are represented 
by second-order accurate centered differences. The difference form of 
the streamwise derivatives is second-order accurate when Crank- 
Nicholson differencing is used and first-order accurate when fully 
implicit differencing is used. The resulting set of coupled, linear 
algebraic equations has a Mock tri-diagonal coefficient matrix with 
4x4 sub-matrices. The equations are solved using a standard tri- 
diagonal matrix inversion method. 
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SECTION ,> 
ANALYSIS 

2.1 COORDINATE SYSTEM 


In this -nalysis, the equations of motion are to be solved using a 
finite difference method. The contours of the duct being studied are 
usually specified in physical space in Cartesian or polar coordinates, 
both designated (z,r). Using physical coordinates for computational 
coordinates, however, is very inconvenient in a finite difference 
method because of difficulties in applying Neumann boundary conditions. 
In addition, the assumptions to be made in deriving a set of equations 
that can be solved by forward marching (see Section 2.2) would, in 
general, be questionable. 

Therefore, an orthogonal, body-fitted, curvilinear system is used 
for the computational U,n) coordinates. The flow domain in the 
physical plane is mapped into a rectangular domain in the computational 
plane, as shown in Figure 2-1. 

Since the governing equations will be written in a general form, 
any orthogonal, body-fitted coordinate system can be used. For the 
results presented in Sections 4.4 and 4.5, an existing computer code 
that employs the coordinate generation method of References 2 and 3 
was used. In this method, the coordinates are generated by solving 
for the planar potential flow through the duct. The computational 
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coordinates are chosen as the velocity potential t and stream 
function n obtained from the potential flow solution. 

The computational and physical coordinates are thus related through 
Laplace's equation, 



a 2 

a n + a n 

? ? 

a z ar 


( 2 . 1 ) 


( 2 . 2 ) 


These equations could be solved directly to get t and n as func- 
tions of z and r. The method used in References 2 and 3, however, 
is to solve the inverse problem, getting z and r as functions of 
t and n* This is accomplished by a conformal mapping technique 
using the Schwartz-Christoffel transformation. 

In order to use this potential flow solution as a coordinate sys- 
tem, the metric scale coefficients hp and h^ must be de- 
termined. From the definition of velocity potential and stream func- 
tion, it can be shown that differential distances along a streamline 
and a potential line are given by 



and 





dn 


respectively, where V p is the potential flow velocity. The dif- 
ferential distance along a general line in this coordinate system is 
then given by 





The governing equations for steady, viscous, compressible, laminar 
or turbulent flow in a two-dimensional orthogonal coordinate system 
are derived in Appendix A as equations ( A. 29 )— ( A. 32 ) . They are repea- 


ted here. 
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ORIGINAL PACE IS 
OF POOR QUALITY 


[( h 2 h 3 Pii) t ♦ (h 1 h 3 pv) n ] . 0 

STREAMWISE MOMENTUM 


(2.4) 
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* ( h l h 3 'E 22 ),] + *E U * ^33 - 'i < 2 ' 6 > 


ENERGY 


F PUCvT 4 + k 0,CvT '' ‘ " W7 p [ (h 2 h 3 u) E * (h l h 3 v) n ] 


- Kfifq [(WEj)* * ( h l h 3lE 2 )n] * *E (2 ’ 7 > 

In these equations, h^, h^, and are the metric scale coeffi- 
cients for the orthogonal curvilinear coordinate system; c and n are 
the coordinate directions; u and v are the velocities in the C and 
n directions, respectively; p, p, and T are the static density, 
pressure, and temperature, respectively; and c y is the specific 
heat at constant volume. The subscripts t and n denote partial 
differentiation. 
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The effective shear stresses are: 

h. 


ORIGINAL PA GC fs 

OF POOR QUALITY 
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r E 22 ■ ^(b 2 \ + 1i^ u )-| *E ' • v 


T£ 33 “ Z “ £ \ 1 h% U * ^ v y- 4 “E » • V 


t E 


12 


, f h 2 /_v_\ + h l/u. 

^2\ h l/n. 




r ( 2 - 8 ) 


J 


" » -■R^0 h 2 h 2 u ' 4 * < h l h 3 v >n] 


For laminar flow, p^. is singly the molecular viscosity p^. For 
turbulent flow, it is an effective viscosity given by the sum of the 
molecular ana turbulent viscosities. 


u t 


+ 


U 1 


where p t » pe, and c, the eddy viscosity, comes from some appro- 
priate turbulence model. 

The effective heat fluxes are: 


q E 1 “ ~ ^7 k E T < 
q E 2 ‘ “ 'RJ k E T n 


(2.9) 
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where, analogous to viscosity, k^ for laminar flow is the molecular 
coefficient of thermal conductivity k L , and for turbulent flow it 
is given by 


k + k. 


where ky is the turbulent thermal conductivity. Here ky is rela- 
ted to uy by 


k 


T 


Vt 


where c p is the specific heat at constant presrure and Fry is 
the turbulent Prandtl number, which is assumed known. 

The effective dissipation 9^, when written out in full, is 
given by: 



Along with the governing differential equations, the perfect gas 
equation of state is used to relate pressure, density, and temperature. 
Thus, 
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p - pRT (2.11) 

where R is the gi'S constant. 

The preceding equations form a set of elliptic, nonlinear, coupled 
partial differential equations to be solved for p, u, v, and T. 

Such a solution is referred to as a complete Navier-Stokes solution 
and can be accomplished using time-dependent or relaxation numerical 
techniques. However, large amounts of computer time and storage are 
required. 

Fortunately, for many flows a complete Navier-Stokes solution is 
not needed. If the flow conditions and geometry are such that elliptic 
effects in one of the coordinate directions can either be ignored or 
assumed known a priori, a marching solution procedure can be used. 
Starting with known profiles at some initial station, the flowfield 
can be computed by integrating the equations one station at a time ... 
the marching direction. The flow field can be computed in one pass, 
and large amounts of computer storage are not needed. 

Therefore, it is assumed that a primary, or streamwise 5 flow di- 
rection can be identified (in this case, the t direction). An 
ore f;r-of -magnitude analysis is used to simplify the governing equa- 
tions (see Appendix B). Viscous and thermal diffusion in the stream- 
wise direction are neglected. Physically, this assumption holds best 
for flows with relatively high Reynolds numbers in which there are no 
abrupt changes in geometry. Separated flow regions cannot be analyzed 
by this technique. 

The above assumptions remove some of the elliptic character of the 
governing equations, but not all. As shown by several autnors (P.efs. 
8-11), the pressure field is also elliptic in nature for subsonic flow. 
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Upstream influence can be important, sometimes very important. Remov- 
ing the elliptic character of the pressure field could result in phys- 
ically unrealistic solutions when using a marching procedure. 

In conventional boundary layer theory, the streamwise pressure 
gradient is assumed to be imposed on the viscous boundary layer by the 
inviscid core flow, and the cross-flow pressure gradient is neglected. 
The imposed pressure gradient usually comes from either an elliptic 
inviscid analysis or experimental data, A similar procedure is used 
in this analysis, but without neglecting the cross-flow pressure 
gradient. 

First, the streamwise pressure gradient p^ is uncoupled from 
the cross-flow pressure gradient p r( as first suggested by Patankar 
and Spalding (Ref. 9). No assumptions are made regarding p^, and 
it is treated as an unknown in the marching solution algorithm. The 
streamwise pressure gradient p^ is written as: 

- P c U*n) ♦ p^U) (2.12) 

where PU,n) is a known two-dimensional pressure field that is "iro- 
posed" onto the flow and p'U) is a correction to P chat is com- 
puted as part of the matching solution to account for viscous blockage 
effects. Note that p'U) has been written as a function of 4 
only, that is, it is assumed constant over any given cross-section of 
the duct being analyzed. The h v 4»n) term can come from any avail- 
able source. Normally, it would be computed using an elliptic invis- 
c id analysis. 

Note that this procedure introduces an inconsistency in the treat- 
ment of the pressure in the two momentum equations. As pointed out in 
Reference 9, this inconsistency is part of the price that must be paid 
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in deriving a set of equations that can be solved by a forward march- 
ing technique. 

With the above assumptions, and the equation of state, equations 
(2.4) -(2.7) can be written in dimensionless form as: 


CONTINUITY 


1 


(w c + sU t> + 


r 2 [‘ h 2 h 3> t ’ u * < h l h 3 ) n e {l ■ ° 


(2.13) 


STREAMWISE MOMENTUM 


r- P uu r + — pvu_ - ~~ pv + nr puv = 


'C h 2 "’“r, h } h 2 


X. + % 


h l h 2 


- -i- <p r * p; ) - 1 


' r c V Re r 


1 / \ % \ 
( h l h 3 T E 12 )n Re r \ h i h 2 E 12 h 1 h 2 E 22 ^3 E 


33 


(2.1*0 


CROSS-FLOW MOMENTUM 


i- puv t ♦ ^ pyv^ + hjh 2 


"1 

n 2 

puv - j-c- p u = - 


EJ R<oT) n + 


( h 2 h 3 T E 12 )c * ( h l h 3 T E 22 ) n] + ^ ^ 


33, 


(2.15) 


ENERGY 


A P<*..T t + Tu »% T . 


- - E^ ’ RT [< h 2 h 3 u >c * < h lV>n] 


Re r Pr r h^h 2 h 2 ( h l h 3 q E 2 ^ n *£ 


(3.16) 
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In these equations, length has been nondimensionalized by L r , 
velocity by U r , density by p r , temperature by T r# viscosity 

p 

by u r , thermal conductivity by k r , pressure by p r U r , and specific 
2 

heat by U r /T r . The reference Reynolds and Prandtl numbers, then, are 
given by 


and 


Pr 


r 



U 2 u 

r M r 

o: 

r r 


When nondimensionalized, the shear stresses (except for the tr in 

t 12 

the cross-flow momentum equation), the heat flux, and the dissipation 

are still given by equations (2.8M2.10) . In the t* term, the 

t 12 

higher-order (v/h^) half has been dropped, since leaving it in would 
make the equation elliptic. As shown in Appendix B, many of the shear 
stress terms in equations (2.14) and (2.15) are actually higher-order 
terms that could reasonably be eliminated. They have been left in for 
completeness, however, since they can be handled within the framework 
of a numerical marching solution. In the computer code used to solve 
these equations, the user can control which viscous terms are included 


in the solution. 
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2.2.2 Boundary Conditions 


Boundary conditions are needed for the dependent variables p, u, 
v, and T at a symmetry line and at a wall. At a symmetry line, the 
boundary conditions are 



At a wall, the conditions normally used are 





(2.17) 


(2.18) 


The density boundary condition is found by writing the cross-flow mo- 
mentum equation at the wall, using the equation of state to rewrite 
ap/»n in terms of p and T. The resulting equation, when written 
in difference form, can be solved for p w in terms of the other 
dependent variables. 

The above velocity boundary conditions are general to allow for 
slip or a moving wall, and to allow for bleed or blowing. For no slip 
and no bleed, u w * v w = 0. The temperature gradient boundary con- 
dition is also general. For an adiabatic wall, the temperature gra- 
dient normal to the wall would be set equal to zero. Details on the 
implementation of these boundary conditions in the finite-difference 
equations are presented in Section 3.3.3. 
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In turbulent flow, when equations (2.18) are used as boundary con- 
ditions at a wall, the finite-difference grid must be very dense near 
the wall to resolve the steep velocity gradients expected there. This 
can be avoided, however, by making use of the experimentally observed 
properties of turbulent boundary layers (Ref. 38). A typical unsepa- 
rated turbulent boundary layer can be divided into an inner region, 
where the total shear stress is essentially constant, and an outer 
region, where the total shear stress decreases with distance from the 
wall. The inner region can be 'urther divided into a laminar sublayer 
where the laminar shear stress oc. ,T ''..ates, a transition or buffer region 
where both laminar and turbulent shear stress are important, and a 
fully turbulent region where turbulent shear stress dominates. 

In the fully turbulent part of the inner region, a universal 
velocity distribution law exists, given by Reference 38 as 

u* - ! In y + + B (2.19) 

K 

where * is the von Karman constant (about 0.4) and B is a con- 
stant between 4.9 and b.b. The variables u + and y + are de- 
fined by 


T 



where y is the distance from the wall, v w is the kinematic viscos- 
ity at the wall, and u t is the friction velocity defined by 
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Equation (2.19) is often referred to as the logarithmic law-of-the- 
wall. It is valid for y + between about 35 and 350 (Ref. 28). The 
upper limit on y + actually depends on the streamwise pressure gra- 
dient, but 350 is a conservative value. If it is assumed that the 
first grid point away from the wall is in the law-of-the-wall region, 
a slip velocity (or wall function) boundary condition can be derived. 
The laminar sublayer is not resolved, thus fewer grid points are needed 
in the n direction. 

To use the law-of-the-wall as a boundary condition, equation (2.19) 
is differentiated with respect to n, giving 


au 



1 

y 


( 2 . 20 ) 


where the + sign is used for a lower wall and the - sign is used 
for an upper wall. This equation is the one actually used to specify 
the boundary condition on u. 

It should be noted that equation (2.19) was developed for incom- 
pressible flow. Very little error is introduced, however, by using it 
for compressible subsonic flow, as will now be shown. Several methods 
have been proposed to extend equation (2.19) to compressible flow 
(Ref. 29). Maise and McDonald (Ref. 30) suggested using equation 
(2.19) with u replaced by a generalized velocity u'. For adiabatic 
flow they used 



u 


( 2 . 21 ) 
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where 


a 




K 


1 


ZT7 

1 "e 


Here \i a and M Q are the velocity and Mach number, respectively, 
at the edge of the boundary layer. For M e « 0 (incompressible 
flow), equation (2.21) gives u' « u. As M g increases, the value 
of u' increases for a given value of u. But, for u/u g ■ 0.4, for 
example, the value of u'/u e is only 0.4018 at M g « 1, which is 
less than a one percent increase. Since the present analysis is for 
subsonic flow, equation (2.19) can be used as is, with no modification 
for compressibility. 

When a wall function boundary condition is used for u, and the 
wall temperature is specified, a wall function boundary condition must 
also be used for temperature since the sublayer region is not resolved. 
(If the temperature gradient normal to the wall is specified, the im- 
plicit assumption that q « q w in the law-of-the- , . j a'' 1 region means 
the temperature gradient boundary condition of equation (2.18) can 
still be used). Following White (Ref. 28), but without assuming 
Pr T « 1, a temperature 1 aw-of-the-wal 1 can be written as 

T + * Pr T - In y + + c, (2.22) 

T K i 


where 


+ 

T 


p w c p u x^ T - v 

% 


c x a 12.8 Pr°* 68 - 7.3 
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Here Pr L « c p y, /k is the laminar Prandtl number and q w is the heat 
flux at the wall. Equation (2.22) is differentiated with respect to n 
to give 





T - T, 


( pr T 7 





(2.23) 


Details on the use of equations (2.20) and (2.23) as boundary condi- 
tions are presented in Section 3.3.4. 


2.2.3 Initial Profi les 


To start the marching procedure, profiles of p, u, v, and T 
must be specified at some initial station. Two options are available 
in the computer program written for this study. A uniform free stream 
can be specified, with boundary layers computed for each wall from 
input boundary layer parameters. Or, complete nonuniform initial pro- 
files of p, u, v, and T can be read in, with the option of computing 
boundary layers at the walls from input boundary layer parameters. In 
either case, the procedure for computing the profiles in the boundary 
layers is the same. 

In laminar flow, the boundary layer thickness 6 and the displace- 
ment thickness 5* are specified. The streamwise velocity in the 
boundary layer is computed by assuming a von Karman-Pohl hausen pro- 
file, that is. 



(2.24) 
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where 

* 

A - 36 - 120 4“ 
0 


Although equation (2.24) was originally developed for incompressible 
flow, it can be used in compressible flow as well (Ref. 31). 

For turbulent flow, the situation is more complicated. Coles 
(Ref. 32) developed an expression for the velocity profile in the 
fully turbulent part of the boundary layer. This was modified by Walz 
(Ref. 31) to include the laminar sublayer and the transitional region. 
Walz' expression is 

u + = = ^1 - i - Ba^ y + - B e" ay + ~ ln(l + y*) + B + ^ u»f* 

(2.25) 



The only symbol in these equations net previously defined is the con- 
stant a, which Walz sets equal to 0.3. 

Before equation (2.25) can be used, the shear stress, viscosity, 

and density at the wall must be found. The boundary layer parameters 

★ 

specified are the thickness 5, displacement thickness 6 , and mo- 
mentum thickness e. The wall shear stress is found from the 
Ludweig-Ti liman skin friction formula, as written by Sasman and Cresci 
(Ref. 33). Thus, 
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T W 

c f m y ~ 2 ■ 0.246 

7 c e u e 



(2.26) 


where T is a reference temperature given by 



i T + 0.22 rv + (0.5 - 0.22 rv) 

^ >o i- h 'o 

e e 


(2.27) 


Tq is the total temperature at the edge of the boundary layer and rp 

® 3 

i s the recovery factor, assumed equal to V’v The parameter v is 

the kinematic viscosity evaluated at T * T. The incompressible shape 

factor H i is computed from 


H * 


'w 




H 1 


M 


(2.28) 


where 


H 


6 . 

6 


The wall temperature is either given, or computed assuming an 
adiabatic wall from 


t aw ■ T e + r F TcJ - (2>29) 

The viscority at the wall is then found from Sutherland's formula, 
equation (2.32). Finally, the wall density is computed from the equa- 
tion of state by assuming ap /ay - 0 across the initial boundary 
layer. 

For both laminar and turbulent flow, the temperature profile in 
the boundary layer is given by the Croce o-Busemann relation (Ref. 28), 
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V 


( J AW “ T W 


u e 




( 2 . 30 ) 


where r f - for laminar flow and ^Fr^ for turbulent flow. 
The density profile is then found by assuming ap /ay - 0 across the 
boundary layer. That is, 

p * (2.31) 
Finally, the cross-flow velocity v is set equal to 0.0. 


2.2.4 Molecular Transport Properties 


In addition to equations (2.13)-(2.16), equations are needed rela- 
ting u L , k L , and c y to temperature. The molecular viscosity 
is determined from Sutherland's formula. 




v 


r 


T r + 110*3 / t V 

t ♦ uo;3 \r^) 


(2.32) 


where T is in K. The specific heat at constant volume is found from 
the empirical formula of Reference 34, which for air with an average 
molecular weight of 28.97 kg/mole can be written as: 

C . 1.44xl0 3 - 3.94xl0 3 r 1/2 - 1.943x10 s T’ 1 ♦ 4.09xl0 7 T -2 - R 


(2.33) 

2 2 

where T is in K and c y and R are in m / sec - K. The 
molecular coefficient of thermal conductivity is given by the 
empirical formula of Reference 35: 
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k 


L 


°QlGlXfi I f. 

°* PQi ' 


a VT 

1 ♦ j x 10” C/T 



(2.34) 


where, for T in K and k in kg-m/sec^ - K, a » 0.002646, b - 245.4, 
and c - 12. 


2.3 STREAMWISE PRESSURE GRADIENT 


2.3.1 Viscous Pressure Correction 


Equations (2.13)-(2.16) form a set of four coupled, nonlinear, 
partial differential equations to be solved for p, u, v, and T. 
However, by rewriting the streamwise pressure gradient in the form of 
equation (2.12), an additional unknown, p^, has been introduced. The 
standard procedure for computing p^ during a matching step is as 
follows (Refs. 9, 20): 

1. Guess p^. 

2. Solve some form of the streamwise momentum equation for the u 
velocity distribution at the new station (after uncoupling and 
linearizing the equation by lagging quantities one step). 

3. Compute the mass flew rate through the duct at the new station 
(which, in general, will be incorrect) using 

f pu dA - m (2.35) 

where A is the duct area and m is the total mass flow rate. 

4. Correct p^ to give the proper mass flow rate. 
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In this analysis, however, the differential continuity equation is 
to be solved simultaneously with the momentum and energy equations. 

The mass flow rate through the duct is automatically conserved during 
a marching step no matter what value of p^ is used, at least within 
the truncation error of the numerical method. Conservation of total 
mass flow rate, therefore, cannot be used to compute p^ after equa- 
tions (2.13 )- (2.16) are simultaneously marched one step. It can still 
be used, however, if p£ is computed before each main marching step. 
Then, once is known, equations ( 2 . 13 )— (2 .16) can be solved simul- 
taneously with p^ treated as a source term. 

In this analysis, in other words, a preliminary marching step is 
taken using a modified form of the procedure outlined above to compute 
p^'. The streamwise momentum equation is uncoupled from the continu- 
ity, cross-flow momentum, and energy equations and linearized by lagg- 
ing certain quantities one step, just as above. The total mass flow 
rate equation is written as 



(2.36) 


The superscript * is used on u to distinguish it from the stream- 
wise velocity computed during the main marching step. Better results 
were obtained using the equation in this form rather than equation 

(2.3 5) . Then, instead of iterating to find the value of p^ that sat- 
isfies total continuity, the streamwise momentum equation and equation 

(2.36) are solved simultaneously as a coupled set to get p^. 

This procedure for computing p^ is not completely rigorous, at 
least not as rigorous as the rest of the solution, since an uncoupled 
form of the streamwise momentum equation must be used. However, p* 
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is only a correction to P, ti.c imposed pressure field, to account for 
viscous blockage. Thus, small percentage errors in p' will give 
much smaller errors in P + p'. In addition, in deriving a set of 
equations that can be solved by forward marching, it was assumed tnat 
the flow is primarily in the i direction. The errors in computing 
the increase (or decrease) in viscous blockage over one marching step 
by using the uncoupled streamwise momentum equation should be small. 


2.3.2 Imposed Pressure Field 


In this analysis, the computational mesh for d gene al case is 
computed using a planar potential flow solution, as described in 
Section 2.1. Thus, for planar two-dimensional flows a potential flow 
velocity can be directly related to the mesh through 

v p-^"^ (2 - 37) 


For axi symmetric flows a separate potential flow solution is required. 
Using the method of Reverence 36, the axisymmetric potential equation. 


( h 2 h 3 \ . / h l h 3 \ n 


(2.38) 


is solved using an iterative alternating-direction implicit method. An 
axisymmetric potential flow velocity is then computed from 



(2.39) 


The "real 14 incompressible inviscid velocity is found from 
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V 


i 


V V 

Vp 


(2.40) 


where 


V 


r 



P AVE 


The subscript AVE here denotes the average value at the initial sta- 
tion. The averaqe incompressible velocity V. comes from 

1 AVE 


y , _ — 

1 AVE p i A l 

where m is the mass flow rate, A^ is the duct area at the initial 
station, and p- is the incompressible, or total, density. For 
compressible flows, the velocity given by equation (2.40) is modified 
using the L ieblein-Stockman compressibility correction, which was de- 
veloped specifically for internal flows (Ref. 37). This compressibi lity 
correction has been shown to yield very good agreement with data over 
a wide range of flow conditions, including high Mach number subsonic 
anc even transonic flow (Refs. 38-42). Good agreement with data was 
obtained in Reference 39 for local Mach numbers as high as 1.3 in small 
renions of supersonic flow. The details of this compressibility cor- 
rection are given in Appendix C. 

The inviscid pressure itself, then, is computed from 



(2.41) 
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where Pq Is the total pressure, a Q Is the *otal speed of sound, 

V c is the compressible velocity, and y is the ratio of specific 
heats. 

It should be noted that, for a given geometry, the potential flow 
needs to be computed only once. The velocity field given by equation 
(2.37) or (2.39) can be saved on a mass storage device. A pressure 
field can then be computed for any initial flow conditions from equa- 
tions (2.40) and (2.41). 

2.4 TURBULENCE MODEL 


The Reynolds stresses and the turbulent heat flux are modeled using 
an eddy viscosity approach to relate these terms to the mean flow (see 
Appendix A). A two-layer model is used. In the outer or wake region, 
the turbulence model of Cebeci and Smith (Ref. 43) is used. In the 
inner or near-wall region, either the model of Cebeci and Smith or of 
McDonald and Camarata (Ref. 44) can be used. The computer code in- 
cludes subroutines for both of these models, which are known as alge- 
braic, or zero-equation, turbulence models. A more complex multi- 
equation turbulence model (i.e., one that solves partial differential 
equations for the turbulence length and/or velocity scales) could be 
incorporated into the code at a later date. However, given the current 
status of multi-equation models and the past success of algeoraic mod- 
els, this probably is not warrented at the present time (Kefs. 29, 4b). 
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2.4.1 Outer Region Model 


In the outer region, the eddy viscosity is given by 


■ aU 6, 

o e k 


(2.42) 


Here, for internal flows, the subscript e implies a value at the 
point of the first maximum streamwise velocity away from a wall. For 
most applications, this is simply the maximum velocity in the duct. 

it 

Also, is the kinematic displacement thickness, 

•y. 




where y is distance from the wall. The parameter o is given by 


1 ♦ fl 

° * °0 T~^ir 


0 


where - 0.0168 and = 0.55. The parameter n is 
Coles' profile parameter, defined as 


n - n^l - exp^-0.243 V^I- 0.2« h)] 


‘ TS Re » k ‘ 1 


He 


e k 

% ■ 


where 
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Here v is the kinematic viscosity and the subscript W denotes a 
wall value. The factor (1 ♦ Dq)/( 1 ♦ fl) is essentially a low 
Reynolds number correction to oq. 

In conventional boundary layer analyses, the eddy viscosity given 
by equation (2.42) ^s often multiplied by an intermittency factor. 

This accounts for the experimentally observed fact that, as the free 
stream is approached, the fraction of time the flow is turbulent de- 
creases. For infer o’ flows, however, this is r.ot necessarily true. 

In fully-developed pipe flow, for example, the eddy viscosity in the 
outer regie : is essentially constant (Ref, 46). Also, in their calcu- 
lation of developing turbulent pipe flow, Richman and Azad (Ref. 47) 
used both a constant outer-region eddy viscosity and a variable eddy 
viscosity that decreased as the centerline was approached. They found 
that the computed mean velocity profiles were insensitive to the outer- 
region eddy viscosity distribution. It was therefore decided not to 
use an intermittency factor in this internal flow analysis. 


2.4.2 Cebeci-Smith Inner Region Model 


In the Cebeci-Smith model, the inner-region eddy viscosity is 
given by 


e 


i 



(2.43) 


where t, the mixing length, is given by 
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» - Ky(l - e~ y/A ) 


(2.44) 


Here « is the von Kanman constant, set equal to 0.40, and y is 
distance from the wall. The parameter A is the well-known Van Driest 
damping parameter. As modified by Cebeci for flows with pressure 
gradient, heat transfer, and mass transfer (Refs. 48-49), A is given by 


A m 



where 


A 
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for flows without bleed (v* - 0), 
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1 - 11.8 
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P 
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2.4.3 McDonald-Camarata Inner Region Model 


The inner-region eddy viscosity in the McDonald-Camarata model is 
still given by equation (2.43), but the form of the mixing length is 
different. They use 

i * 6 tanh 9 (2.45) 

where = 0.09, * * 0.43, and 6 ■ y at the point where 
u = 0.99 u g . The damping factor 3 is given by 

9 = ^ 1/2 q ~ j (2.46) 


where & is the normal probability function, y = 23, and = 8. 
The parameter y + is found from 



where 


u 


T 



Note that in this model, the friction velocity u t is computed 
from local values of density and shear stress, not wall values. 

The boundary between the inner and outer regions is determined by 
requiring that the eddy viscosity be a continuous function of distance 
from the wall. In other words, starting at the wall, e = 

> Eq, then e « e q . 


until 
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For ax i symmetric flows with a centerbody and for nonsymmetric 
planar flows, the eddy viscosity is computed separately for each wall. 
In general this leads to a discontinuity in the two outer-region val- 
ues. In this case, the eddy viscosity is made continuous by linearly 
interpolating between the two values over the outer 50 percent of each 
outer region. If never reaches e 0 for either boundary 
layer, the point of maximum is used. In addition, the interpo- 
lation region is forced to include at least the outer 50 percent of 
each boundary layer. 

An example of a typical eddy viscosity profile is presented in 
Figures 2-2a and b. This particular profile is for developing turbu- 
lent flow in a pipe. Figure 2-2a shows the entire inner region, with 
e- increasing from zero at the wall to the constant outer region 
value cq . In Figure 2-2b, the profile very near the wall is pre- 
sented, showing the damping effect in the viscous sublayer. 
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SECTION 3 

NUMERICAL METHOD 

3.1 BASIC DIFFERENCING PROCEDURE 


3.1.1 Computational Mesh 


The governing equations presented in Section 2 are solved using a 

♦ 

finite-difference method. A computational grid, or mesh, is thus re- 
quired in the U.n) coordinate system. Since the method is to be 
applied to viscous flow problems, grid points should be closely spaced 
near the walls, where the velocity gradients in the n direction are 
expected to be largest. For convenience when applying difference 
formulas, however, a uniformly spaced mesh is desired normal to the 
walls. The desired grid is obtained by employing the following trans- 
formation (Ref. 2): 


n 



(3.1) 


C - \ (1 - 2D)" 1/2 


where 
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and D is given by 


1 



-0 In 


D 

7 


In equation (3.1), Y is the equally spaced, or transformed, coor- 
dinate normal to the walls and n is the packed, or untransformed, 
coordinate. Note that Y is a function of n only. The transforma- 
tion is implemented by using chain rule differentiation in the form 

a _ dY a 
3n dn 17 


where dY/dn can be derived from equation (3.1) as 


dY 1 

^"7 




In 



(3.2) 


For geometries with two walls (e.g., an axisymmetric flow with a cen- 
terbody or a nonsymmetric planar flow), both Y and n in equation 
(3.1) vary from 0 to 1. Then, for uniform spacing in Y, the trans- 
formed coordinate, equation (3.1) will pack the n points close to 
both walls. For cases with only one wall, Y, and thus n, are allowed 
to vary only from 0.5 to 1. The n points are then packed only at 
the outer wall. This mesh is then expanded down to the centerline, 
where Y * n = 0. The degree of grid packing is controlled by the 
parameter D^, which would typically be between 5 and 500. Increas- 
ing packs more points near the walls. 

It should be noted that no transformation is used in the streamwise 
direction, since the locations of high gradients in that direction will 
vary from case to case. However, the marching step size aC is not 
assumed constant. The user can therefore tailor the distribution of 
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grid points in the t direction to the particular case being 
studied. 

The relationship between the U,n) and (c,Y) coordinate systems 
is shown in Figure 3-1. Also shown in the figure are the indices used 
in the finite-difference grid. An "i" is used as the index for the 
t, or marching, direction, with i = 1 at the initial station. A 
"j" is used in the cross-stream direction, with j » 1 at the lower 
wall or symmetry line and j = J at the upper wall. 

3,1.2 Difference < ormulas 


Standard difference formulas can be derived using Taylor series 
expansions. The derivations can be found in any basic reference on 
finite-difference methods (e.g.. Refs. 50-52). Therefore, only the 
resulting formulas are presented here. 

First derivatives in the marching direction are represented by 



(3.3) 


This expression is either first or second order in t, depending on 
the value of the implicit weighting factor w. This will become 
clearer in Section 3.2. 

Derivatives in the Y direction at interior points (2 £ j <_ J-l) 
are all represented by standard second-order centered differencing 
formulas. For first derivatives the formula is 



F 


J+l ~ F J - 1 
“ 2 iY 


(3.4) 


and for second derivatives it is 
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(3.5) 


Mixed second derivatives are simply given by a combination of 
equations (3.3) and (3.4). That is. 


(_& ) .i. 

VKlY/ (tWiJ At 


(4 y F i*i,j 


- *V F i,j> 


(3.6) 


Difference formulas for Y-deri vatives are also needed at boundary 
points (j « 1 and j - J) in order to apply some of the boundary con- 
ditions described in Section 2.2.2. These formulas are presenteo in 
Section 3.4. 

With the above difference formulas, the resulting set of coupled 
algebraic equations will have a block tri-diagonal coefficient matrix 
with 4x4 submatrices. 


3.2 LINEARIZATION PROCEDURE 


When a marching step is taken from station i, where the solution 
is known, to station i+1, where it is unknown, the finite-difference 
equations are written at the intermediate station i+w. Here w is 
the implicit weighting factor, 0 < w < 1. For w « 1 the difference 
formulation is therefore fully implicit, while for w ■ 1/2 it is of 
Crank-Nicholson type. This explains why the difference formula for 
E-derivatives (eq. (3.3)) can be either first or second-order. For 
v: * 1/2 it is a second-order centered difference formula, while for 


w » 1 It Is a first-order backward difference formula. Different 
values of w can be used in each of the four governing equations. 

When this type of implicit weighting is used on a simple linear model 
equation, stability analyses show the method to be unconditionally 
stable for 1/2 ^ w < 1, but only conditionally stable for 0 < w < 

1/2 (Kef. 52). For tie complex system of nonlinear equations consid- 
ered in th's analysis, however, these stability conditions can only be 
used as guidelines (Ref. 51). The results to be p .sented in Section 
4, therefore, were computed using w - 1 throughout because of its 
greater margin of stability. 

When the governing equations are written at i+w, the unknown 

terms are, in general, nonlinear functions of the unknown dependent 

variables p, u, v, and T. For example, in the streamwise momentum 

9 U d u 

equation the terms pu — and pv -^y appear. These nonlinear terms 
must be linearized in order to solve the final set of algebraic equa- 
tions by matrix inversion. 

One method of linearization would be to simply evaluate the coef- 
ficients at the known station i and only evaluate the derivatives at 
i+w. However, this is only a first-order method, unless iteration is 
applied. In the present analysis, a second-order linearization proce- 
dure is desired, but without requiring iteration. 

This is accomplished by expanding each unknown nonlinear term in a 
Taylor series in c about the known station at This same proce- 

dure has been used in other finite-difference marching analyses (Refs. 
2, 53-54). Let G * G(p,u,v,T) represent a general nonlinear term. 
Expanding in a Taylor series gives 
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G in - G > * ('ll), « * 0 («> 2 

where At - and 

aG + + 

at “ ap at au at av at aT at 

Note that for G.^ to be second-order accurate, J and therefore 
(f), etc., need only be first-order accurate. Note also that a 
general term G will include metric scale coefficients and the grid 
stretching parameter dY/dn. However, these are known a priori and can 
be evaluated at i+w directly. 

In the governing equations, there are six basic forms of nonlinear 
terms, as follows: 



Here f is a function of p, u, v, and T. In the second and third 
basic forms, Q is also a general function of p, u, v, and T. In 
the last three, however, it is identically equal to p, u, v, or T. 

The first three forms, which do rot include ^-derivatives, are 
written at station i+w as 

6 i+w " "S+l + (1 “ w > G i {3 ‘ 9) 

(For all equations in Section 3, terms without an "i" or "j" subscript 
are understood to be at grid location i or j, unless stated 
otherwise) . 

App ying equations (3.7) and (3.8) to the first basic expression, 
simply f, gives 


( 3 . 7 ) 


(3.8) 
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Vi ■ f i * (if), « 


f 


i+l 


f i* 


/ 8f ap . af 9u . af »v . af aT \ Ar 
\a 0 H au as av ic TT as ^ 


Then, using forward difference formulas for the s-derivatives and 
applying equation (3.9) gives 

Vw ■ f t* w (ff lp *lf 4u *'J7 4v+ lf 4T ), < 3 - 10 ) 

where (a f ). * p. +i - p., etc. The nonlinear term f^ is thus 
represented as a linear combination of the dependent variables at the 
unknown i+l station. As an example, using a term from the stream- 
wise momentum equation, let 


f 



CUV 


The subscript n denotes partial differentiation. Writing the 
metrics at i + w directly, and using equation (3.10), 

{'iVt * w [Vi<°in " * VdVi - u i> 

i+w 



+ 


'i u i< v 


i+l 



The second basic expression is f where both f and Q can 
be functions of p, u, v, and T. Equation (3.7) gives 


Now 


(-a),., - 
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M'S) 

Combining, and using equation (3.b), gives 

('S)„- (-S),* (SMSS-SMSS),®*- 

a /^^ + iaau + aQ3v + iqiI\ At 
i 7? \ 3 p at au at av at aT at ^ 


+ f. 


Finally, using forward differences for the t-derivatives and apply- 
ing equations (3.9) and (3.4), 

af ... . af 

, * 4 "tu AU T 

i+w 


(f TyL. * {U ^ ] i 4 w Ap + Su Au * iv AV + TT «)** 

*» * w lu + lv sv lT ), t3-U) 


+ Wt 


The third basic expression is |y (f ~y)* It* linearized and 
differenced form can be written directly using equation (3.11). 

[It (f $)] * «,[lf ‘ Y ») 1 ♦ «(£ “> * & lu * w iv * It 

1 w 

* * i { Y (h 10 + !u Lu * « tV * $ ‘ T ), * w(4 Y f ) i*Y (l7 10 


♦ itQ. au + av 
au av 




(3.12) 


The fourth expression to be linearized is f -|^. This is written 
at i+w as 



i+w 




At 


i 
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Recall that for the last three basic expressions to be linearized, Q 
is identically p, u, v, or T. To evaluate f i+w , equations (3.7)- 
(3.9) are used. That is, 

f 1*w - "Vl * - "> f i 

where 

f i+l * f i * (If), « 


However, to maintain linearity of the difference equations, (-^) must 
not contain terms at i+1. Therefore, backward differencing is used in 
this case. Thus, 



where “ 5 i " *i-l' 


Then 


i + w 


wxfr - 


f i-l> 


where 



Finally, the original expression becomes 



V *< f i 




(3.13) 


The fifth expression to be linearized is (f Its linear- 

ized and differenced form can be easily written using equation (3.13). 


[I* ( f f )],./ [v «»<'i 



+ 



at 


{y [f, ♦ 


W»(f, - f w l] 


(3.14) 
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Finally, the sixth expression is ~ (f -$)• This is written, 
using equations (3.13) and (3.8), as 

, - 


(« V Q) 


Y^'i+w 


L U Vw 


i+1 


- Pw 


a* 


V#u Vw 


i+1 


- u. 


AC 



(3.15) 


where the combination 



(t?L is 9lven b¥ 

1 + W 



etc. 

Every term in the governing equations was linearized and differ- 
enced using equations (3.10)— (3- lb) , with two exceptions. The first is 
the effective viscous dissipation in the energy equation, given by 
equation (2.10). While this term could also be linearized to second- 
order in c, the resulting difference form would be quite complex. 

It was therefore decided during the initial code development to use a 
first-order linearization, evaluating the viscous dissipation at the 
known station i. For the cases computed to date, this has proved to 
be sufficient. 

The second exception is the turbulent viscosity py (and there- 
fore the turbulent thermal conductivity ky, also). The turbulent 
viscosity is too complex a function of the dependent variables to lin- 
earize using the formal procedure presented in this section. It is 
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therefore computed using the turbulence model equations of Section 2.4 

with the dependent variables evaluated at the known station. The 1 am- 

inar viscosity m l , however, is treated formally as a function of 

temperature during the linearization procedure. 

The difference equations resulting from the application of this 

1 inearization procedure are quite long, even when the difference oper- 
2 

ators 6y and are used. Therefore, they are not written here, 
but instead are presented in Appendix D. 

3.3 BOUNDARY CONDITIONS 


3.3.1 Derivatives at Boundaries 


In order to solve the difference equations, difference forms of 
the boundary conditions presented in Section 2.2.2 are needed. In par- 
ticular, a difference formulation for where Q is one of the de- 
pendent variables, is needed at j * 1 and j ■ J. A three-point, 
second-order, one-sided difference formula can be derived using Taylor 
series. Thus, 



+3Q ± 
, jL 


4 v 



(3.16) 


where the top sign is used for the j « 1 boundary and the bottom sign 
is used for the j ■ J boundary. Solving for Q- gives 

J 

"j -j[ T(24V )(^)/ 4 ^l ' «j«] 


(3.17) 
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With this formula, the tri-di agonal ity of the coefficient matrix is 
retained. 

When the cross-flow momentum equation is evaluated at a wall to 
derive the density boundary condition, a one-sided difference formula 

3 2 V 

is also needed for — *. Unfortunately a second-order one-sided d if— 
ference formula would require more than three points and thus destroy 
the tri-di agonal Uy of the coefficient matrix. Thus, a first-order 
formula was employed, given by 



(3.18) 


where, again, j * 1 or J. Note that this first-order formula at the 
wall is the same as the second-order formula one point away from the 
wal 1. 


3.3.2 Symmetry Line 


At a symmetry line, the boundary conditions are given by equations 
(2.17). Using equation (3.17) with j ■ 1, these are expressed in dif- 
ference form as 


P 1 

U 1 


3 (4p 2 
3 ( 4u 2 
3 (4T 2 









v 


1 


- 0 


J 


(3.19) 


where the subscripts refer to j-locations in the grid. 
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3.3.3 Solid Surfaces 
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At a solid surface, the usual boundary conditions for u, v, and 
T are given by equations (2.18). The conditions for u and v in 
difference form are simply 

u j - u w (3.20) 

v 0 - v w (3 - a) 

where j * 1 or J, and u w and v w are known functions of t. For 
no slip and no bleed, u^ = v w - 0. The condition for specified 
wall temperature is 

T j - T w < 3 ' 22a) 


For a specified temperature gradient at the wall, equation (3.17) 
gives 


T . 
J 



(3.22b) 


For an adiabatic wall, of course. 



= 0 . 


The density at a solid surface is found by evaluating the cross- 
flow momentum equation at the surface. The difference equation at the 
wall is the same as the difference equation at interior points, and is 
presented in Appendix D as equation (D.3). At the wall, however, the 

o 

6y and 6y operators in equation (D.3) are given by equations (3.1b) 
and (3.18). When the resulting algebraic equation is solved for p^. 


the result has the form 
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>j ‘ C W (S’j*! * C W 3 P j*2 + \ U j + C W 5 u j*l + C W & u j*2 * C W ? V 


4 C„ v JJL , + C„ v^ 0 ♦ C., 1 . ♦ C„ T, a1 ♦ C Ll T^o ♦ ^(3.23) 


W 8 >1 ^>2 > n j Wn 3*1 L W 19 'j*2 


'10 


11 


’12 


where, again, j ■ 1 or J. Details on the derivation of this equa- 
tion are presented in Appendix E. 

The boundary conditions on p, u, v, and T that are normally 
used at a solid surface are thus given by equations (3.23), (3.20), 
(3.21), and either (3.22a) or (3.22b). 


3.3.4 Wall Functions 


In turbulent flow, if it is assumed that the streamwise velocity 
near the wall obeys the logarithmic law-of-the-wa II, a slip velocity 
(or wall function) boundary condition can be derived. The velocity 
gradient normal to the wall in ...e law-of-the-wa 11 region is given by 
equation (2.20). Applying this equation, in difference form, at the 
first point away from the wall gives 


j*2 


- L 





1 



Solving for Uj gives 


u j - u j*2 * ( * 4V ) 




(3.24) 


where j = 1 or J, and y^ is the distance from the appropriate 

wall. The friction velocity u is computed from the results at the 

T j 

previous station by a secant iteration procedure. 
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If a wall function boundary condition is used for u, it must 
also be used for T if the wall temperature T w is specified. 
Applying equation (2.23), in difference form, at the first point away 
from the wall gives 


T j*2 T j / dY \ 

liV l dr, A 


T j*i - T W 


j*l 2 j*l <yj*i(Pr T i In y * 41 ♦ Cj) 


Solving for T. gives 

J 


T .w - 


- T j*2 (2iY) (w^) n — ( Pr J T *i y * tl ♦ Cl ) 13-2 ’ 


where again j = 1 or J, yj*i is the distance from the appropriate 

♦ V 


wall, and y = . 

V W 

The other conditions are the same as those given in Section 3.3.3. 
Thus, when wall functions are used at a solid surface, the boundary 
conditions on p, u, v, and T are given by equations (3.23), (3.24), 
(3.21), and either (3.25) or (3.22b). 


3.4 TRI-DIAGONAL MATRIX INVERSION ALGORITHM 

When the boundary conditions are applied to equations (D.5)-(D.8), 
the resulting equations can be written as 


B 2 X 2 + C 2 X 3 “ S 2 


Vj-l * Yj - E J S J*1 ■ U < j < J - 2) t (3.26) 


A J-l X J-2 + B J-l ll J-l ' S J-1 
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mm mm 

where the A's, B's, and C's are known 4x4 matrices whose elements 
are the coefficients of equations (D.5)-(D.8), the Vs are 4-element 
vectors made up of the source terms in equations (D.5)-(D.8), 



and the X's are the 4-element solution vectors, 

p 

u 
v 


L T j 


The A's, B's, C's, and S‘s are presented in Appendix F. Equations 
(3.26) can be written in matrix form as 



(3.27) 
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The equations thus have a block tri-diagonal coefficient matrix whose 
elements are 4x4 sub-matrices, and can therefore be solved using a 
standard tri-diagonal inversion algorithm (Refs. 51, 55). The algorithm 
for a block tri-diagonal matrix is the same as for a scalar matrix, but 
with scalar divisions changed to multiplications by matrix inverses. 

With the indexing as used in equation (3.26) and (3.27), the 
algorithm is as follows. First let 


F x . B" 1 ^ 


(3.28) 




(3.29) 


Then, for j increasing from 3 to J-l, compute and store 


Vi * ( v 

l i-i ‘ ( 5 J * Vj -*)' 1 ( l i - "‘A-a) 


(3.30) 

(3.31) 


The X's themselves are computed by back substitution. First 


X J-1 * 6 J-2 


(3.32) 


Then, for j decreasing from J-2 to 2, 


X. * G . . - F • -iX * i . 
J J-l J-l J + 1 


(3.33) 


Finally, Xj and Xj are found by using the boundary conditions. 
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3.5 VISCOUS PRESSURE CORRECTION 
3.5.1 Streamwise Momentum Equation 

The viscous pressure correction *s computed by solving the stream- 
wise momentum equation, equation (2.14), and the total mass flow rate 
equation, equation (2.36), during a preliminary marching step, as 
described in Section 2.3.1. In order to do this, the streamwise momen- 
tum equation must be uncoupled from the differential continuity, cross- 
flow momentum, and energy equations. The uncoupling is done by treat- 
ing u and p^ as the only unknowns, evaluating o, v, and T at the 
previous station. The equation is differenced using the difference for- 
mulas of Section 3.1.2, and the unknown terms are written at i+w 
using equation (3.9). The resulting difference equation is fairly 
long. It is presented in Appendix 6 as equation (6.1). 

After collecting terms, equation (G.l) can be written as 

b 2 u* * c 2 u* + d 2 p‘ - 

4 Vj * C jVl * d j p ( ■ (3 < j < j - 2) (3.34) 

a J-l u J-2 * b J-l u J-l * d J-l p t " S X J _ 1 

* 

The "i+l" and "i + w" subscripts have been omitted from che unknowns u 
and p^, respectively. Here the a's, b's, c's, and d*s are the 
known coefficients and the s£s the known source terms that come from 
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collecting terms In equation (G.l). They are written out In full In 
Appendix G. The superscript * Is used on u* to distinguish it 
from the streamwise velocity computed during the main march 4 '’ '] step. 


3.5.2 Total Mass Flow Rate Equation 


The total mass flow rate equation, equation (2.36), is written at 
station i+1 as 



P i*l 


U»n) ' 

RT 


Pi 


i+1 


(O u * 


i+1 dA i+l 


m 


(3.35) 


Note that the temperature has been lagged one step, just as in the 
streamwise momentum equation. The differential area is given by 


dA 


h 2 h 3 

dY/dn 


(2i) u dY 


(3.36) 


where u r 0 for planar flow and 1 for ax i symmetric flow. The pres- 
sure correction p' is related to the gradient p p * by 




At 


Substituting the above expressio* into equation (3.35) gives 



p . . . + p ! + p ' At 

11 1 c i+w 

rt: 


u i+l ^i+l 


m 


(3.37) 
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The nonlinear term p' u, +1 has thus been introduced. This term 

*i*w 1 1 

is linearized following the procedure described in Section 3.2, but for 

* 


a general function of p^ and u , instead of p, u, v, and T. The 
result is 


■ ■w* * ^ ( p ^„ - u 0 


:+w 


l-i+W 




Note, in this equation, that the streamwise velocities at the known 
station i are those from the previous completely-coupled main march- 
ing step. Substituting into equation (3.37) and using equation (3.3b), 



(3.38) 


The unknowns in this equation are u. +1 and p' 

1 5 i+w 

When equation (3.38) is integrated numerically using Simpson's 

4» 

rule, and the boundary conditions on u are applied, the resulting 
equation can be written as 


e 2 u 2 ' e 3 u 3 * ' ' 


' * e J-2 u J-2 ’ e J-l u J-l * % ’hi < 3 - 39 > 


where again the "i+1" and "i+w M subscripts have been omitted. Details 
on the derivation of this equation are presented in Appendix (*. 
Equations (3,34) and (3.39) represent a system of J-l equations in 
J-l unknowns (J-2 u ■' s a, id p'). 
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In matrix form, equations (3.34) and (3.39) are written as 



(3.40) 


The system thus has a scalar coefficient matrix that is "a’most tri- 
diagonal". An efficient solution algorithm can be derived using 
Gaussian elimination, in the same way the standard tri-diagonal in- 
version algorithm is derived. With the indexing as used in equation 
(3.40), the resulting algorithm is as follows. First let 


= c 2 /b 2 

(3.41) 

! s X 2 /b 2 

(3.42) 

■ V b 2 

(3.43) 


Then, for j increasing from 3 to J-l, compute and store 
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F j-1 ■ c j /(b j - 

(3.44) 

G j-1 ' ( S *j - 4 0 G j-2)/( b j - a j F j-z) 

(3.45) 

H j-1 ‘ ( d 0 ' a 0 H j-2 ) /( b j ' 4 /o-0 

(3.46) 

6 j * *j ' F j-2 e j-l 

(3.47) 

f * f - H j— 2 e j— 1 

(3.48) 

M M j-2 j-1 

(3.49) 

In equations (3.47)-(3.49), and in equations (3.50) and 

(3.51) below. 

the “equals" sign is used in the FORTRAN sense to mean “ 

is replaced 

by”. Next, compute 


f = f “ H j-2 e J-l 

(3.50) 

S M = S M “ G J-2 e J-l 

(3.51) 

Then 


H - S M /f 

(3.52) 

U J-1 “ G J-2 H J-2 P { 

(3.53) 

Continuing the back substitution for j decreasing frorr 

i J-2 to 2, 


u j g j-i- Vivi - h m p ; (3 - 54) 


* * 

Finally, and u^ are found by using the boundary conditions. Note 
that, for this application, the algorithm can be stopped after equation 
(3.52) since only p^ is needed for the main marching step. 
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The viscous pressure gradient correction computed by this 
procedure tends to oscillate for the first few marching steps. This is 
probably due to starting the calculation with the physically unrealistic 
value of p£ » 0. Under many flow conditions, these oscillations tend 
to damp out as the calculation proceeds. In some cases, however, they 
can become severe enough to stop the calculation. To prevent this from 
happening, the value of p^ was " underrelaxed" using 

<PJ>1*W- <!- Ww** Wl- l3 - 55) 

where (pp*^ is the value computed using equation (3.52), and f R is 
the relaxation factor. It was found that with f R = 0.1 the oscilla- 
tions in p^ damped out within the first several marching steps. In 
addition, the computed results with and without t>nde»relaxation were 
essentially the same for cases in which the oscillations would have 
damped out naturally. The value of f R was therefore set equal to 0.1 
for all the test cases presented in Section 4. 


FJgun y L * R«Utkmsh!p tatwain untransformtd (C,n>«od transform* (f Y) omputottonil 
coordinates, 
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SECTION 4 
RESULTS 

To validate the marching analysis described in Sections 2 and 3, 
several test cases were run. In this section, computed results are 
presented and discussed for five cases: (1) developing laminar flow in 
a circular pipe; (2) laminar flow in a two-dimensional converging 
channel (Jeffery-Ha.,iel flow); (3) developing turbulent flow in a cir- 
cular pipe; (4) turbulent flow in a two-dimensional S-duct; and (5) 
turbulent flow in a typical subsonic diffuser for a supersonic inlet. 
For all cases, the results are compared with experimental data and/or 
exact solutions. 

4.1 LAMINAR DEVELOPING PIPE FLOW 


A basic test case for a subsonic viscous internal flow analysis 
is laminar developing flow in a circular pipe. Experimental data suit- 
able for comparison with analysis are plentiful (Refs. 56-58). Other 
numerical results are also available (Ref. 55). In addition, the com- 
puted results should approach the exact Poiseuille solution far down- 
stream of the pipe entrance. 
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4.1.1 Gr^ and Initial Conditions 

The obvious orthogonal coordinate system for a straight pipe is a 
cylindrical coordinate system. Therefore the metric scale coefficients 
used were simply 


For this coordinate system, the computational U,n) and physical 
(z,r) coordinates coincide so that C « z ana n - r. The configu- 
ration is illustrated in Figure 4-1. 

The imposed pressure PU,n) was set equal to a constant every- 
where in the duct, since for this case there are no elliptic effects 
due to geometry. The streamwise pressure gradient is then completely 
determined by the viscous correction p£(0. This case is thus an im- 
portant test of the method used to compute p£ during the marching so- 
lution. 

The entrance length L g for developing pipe flow (i.e., the 
distance from the entrance at which the flow can be considered fully 
developed) is given by (Ref. 59) 

L e 

■p — ■ 0.08 Re^ 

where 0 is the the pipe diameter and Re^ is the Reynolds number 
based on D and Uq, the average velocity in the pipe. The pipe diam- 
eter and flow conditions were chosen to give a Reynolds number low 
enough to allow the computation of the entire entrance length within a 
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reasonable number of pipe diameters. A pipe radius of 0.061 cm (0.002 
ft) was used as the reference length L^. The reference velocity U r 
was 6.10 m/sec (20 ft/ sec). This corresponds to the centerline veloc- 
ity for fully-developed flow, which is twice the average velocity in 
the pipe. Standard atmosphere values of 288 K (519* R) and 1.2246 

O 0 

kg/m (0.07645 lb m /ft ) were used for the reference temperature T r 
and density o r , respectively. The reference Reynolds number Re f was 
thus 254.151, which corresponds to an entrance length of 40.66 radii. 

The calculation could not be started exactly at the pipe entrance 
because the u and v velocity profiles would be singular at the 
wall. It was therefore started slightly downstream of the actual en- 
trance, at z/R = 0.254. The initial streamwise velocity, nondimen- 
sionalized by the average velocity Ug, is presented in figure 4-2. 

This profile was given by the tabulated finite-difference results of 
Hornbeck (Ref. 55), with weighted average quadratic interpolation used 
to get values between the tabulated points. The interpolation resulted 
in a slight kink in the profile near the edge of the boundary layer. 
This kink was quickly eliminated, however, through viscous effects and 
had no appreciable effect on the computed results. The initial cross- 
flow velocity v was set equal to zero. It should be noted that this 
is not physically realistic. In fact, the v-velocity should be 
largest near the pipe entrance. However, it should still be small 
compared to the u-velocity. In addition, it was felt that setting v 
equal to zero was representative of the way this analysis would be 
used in a practical situation, where the initial v-velocity profile is 
usually unknown. Finally, since the flow in this case is essentially 
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incompressible, the initial dimensionless density and temperature were 
both set equal to one. 

No-slip and no-bleed boundary conditions were used for u and v 
at the wall. An adiabatic wall boundary condition was used for the 
temperature. Fifty-one grid points were used between the centerline 
and the wall. The points were lightly packed near the wall using the 
transformation given by equation (3.1) with * 5. The marching step 
size at was 0.05. This grid is plotted in Figure 4-3. Since the 
same grid was used through the entire pipe, only a short section near 
the entrance is shown. To reach the end of the estimated entrance 
length at t * 40.66, 809 marching steps were required. The calcula- 
tion took 3.5 minutes of CPU time on an IBM 370/3033. It should be 
stated, however, that a mesh size study was not done for this case. 
Based on the results to be shown in the following section, it is felt 
that as could be increased considerably, especially downstream of the 
region very near the entrance. The number of transverse grid points 
could probably also be reduced. 

4.1.2 Results 

The computed axial velocity profiles are presented in Figure 4-4. 
Profiles are plotted every ten marching stations, with the abscissa 
displaced by 1/4 of the interval width for each successive profile. 

The boundary layer thickness increases rapidly near the entrance and 
soon reaches the duct centerline. It takes much longer, however, for 
the centerline velocity to reach its fully-developed value. 
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The computed cross-flow velocity profiles are presented in Figure 
4-5. Although the initial v-velocity was set equal to zero, the anal- 
ysis quickly generates a realistic profile. The computed v-velocity 
is negative, with the flow moving away from the wall, and its magnitude 
is a maximum near the edge of the boundary layer. As the boundary 
layer thickness increases, the point of maximum v-velocity moves out 
from the wall. Further downstream, where the flow in the duct is com- 
pletely viscous, the v-velocity is essentially zero. 

The computed u-velocities are compared with the numerical results 
of Hornbeck (Ref. 55) in Figure 4-6. Note that the abscissa in Figure 
4-6 is z/(R Re R ), where Re R = pu Q R/u L . This removes the Reynolds 
number dependence from the problem, allowing the present results to be 
compared directly with other results or data. Hornbeck solved the 
incompressible streamwise momentum equation, minus the diffusion terms 
involving v, along with the total mass flow equation, for the stream- 
wise velocity distribution and the pressure. The agreement between 
the present results and those of Hornbeck is excellent. The velocity 
near the centerline increases continuously until the flow becomes 
fully-developed. Near the wall the velocity decreases continuously. 
In-between (at r/R = 0.6 and 0.7) the velocity first increases to 
compensate for the increasing boundary layer blockage. Then, as the 
boundary layer thickens further, including these points, the velocity 
decreases gradually to its fully-developed value. 

The present results are compared with the experimental data of 
Pfenninger (Ref. 56), Reshotko (Ref. 57), and Nikuradse (Ref. 58) in 
Figures 4-7 to 4-9, respecti vely. In these figures the streamwise 
velocity ra f io u/uq is plotted as a function of z/(R Re R ) at six 


70 


radial locations. Note that Pfenninger's data are all near the duct 
entrance and thus an expanded abscissa is used in Figure 4-7. The 
agreement with Pfenninger's data is excellent. The agreement with 
Reshotko's data, presented in Figure 4-8, is also excellent, except at 
r/R ■ 0.8 and 0.9. However, the accuracy of the data at these radial 
locations, as plotted in the figure, is questionable. The data in Ref- 
erence 57 are presented as velocity profile plots at various a.“ial sta- 
tions. The disagreement in Figure 4-8 can be explained by assuming that 
the data near the wall were actually taken at slightly lower values of 
r/R than 0.8 and 0.9, either intentionally or because of a systematic 
error in determining the probe location near the wall. (The radial 
locations are not explicitly stated in Reference 57, and the plots can 
only be read to an accuracy of about 0.02 in r/R). As an indication 
that this may be the case, note that the experimental velocity near 
the wall does not approach the Poiseuille profile value. (The computed 
velocity does, as will be shown in a later figure). The agreement be- 
tween the analysis and Nikuradse's data is excellent in the downstream 
half of the pipe, but not as good near the entrance. However, these 
data may also be questioned, since they do not agree with the results 
of Reference 55, or the data of References 51 and 52. The present re- 
sults are also compared with the data of References 56-58, in the form 
of velocity profiles at various axial stations, in Figures 4-10 to 
4-12, respectively. 

In Figure 4-13, the computed axial pressure gradient, in the form 
<Kp'/»Uq) 

of (j( z /R is P 1ot ted as a function of z/(R Re R ). Also shown are 

R 

the results of Hornbeck (Ref. 55), the data of Reshotko (Ref. 57), and 
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the fully-developed Polseuille value. Recall that since the Imposed 
pressure is a constant for this case, the axial pressure gradient is 
determined completely by the viscous correction p^. The starting 
transient near the entrance is a result of starting with the physically 
unrealistic (but, as with the cross-flow velocity, realistic in prac- 
tice) value of - 0. The pressure gradient rapidly reaches the cor- 
rect value, however, and slowly approaches the fully-developed value 
far downstream. The agreement with the results of Reference 55 and the 
data of Reference 57 is excellent, indicating that the method used to 
compute the viscous pressure correction p^ is valid. 

In Figure 4-14 the computed skin friction coefficient is plotted 
as a function of axial distance. The shear stress at the wall was 
found by second-order one-sided numerical differentiation of the com- 
puted axial velocity. Also shown in the figure is the value for fully- 
developed Poiseuille flow. As one would expect, the skin friction 
starts out high, where the boundary layer is thin, decreases rapidly 
as the boundary layer thickens, and then slowly approaches the fully- 
developed Poiseuille value. 

Finally, in Figure 4-15 the computed streamwise velocity profile 
at the last station is compared with the exact Poiseuille profile. 

The agreement is excellent, indicating that a ful ly-develped state has 


been reached. 
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Figure 4-1 - Comparison of computed streamwiss vtloclty with date of Rtshotko for laminar 
dsvtloping pip* flew. 



Figure 4-9. - Comparison of computed streamwlsa vtlocity with data of Nikuradsa for laminar 
drreloc'ng pips flow. 
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figure 4-ia - Comparison of computed velocity profiles with data of Pfennings for laminar developing pipe flow. 



Figure 4-U. * Comparison of computed velocity profiles with data of Reshotko for laminar developing pipe floe. 
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4.2 JEFFERY -HAMEL FLOW 


Another useful test case is laminar incompressible flow in a two- 
dimensional wedge-shaped channel, known as Jeffery-Hamel flow. When 
self-similarity is assumed, an exact solution to the Navier-Stokes 
equations exists for this flow. The streamlines are radial, inter- 
secting at a line source (for diverging flow) or sink (for converging 
flow). The solution for the velocity, derived in terms of Jacobian 
elliptic functions, is given by Millsaps and Pohlhausen (Ret. 60). 

4.2.1 Grid and Initial Conditions 


The configuration and coordinate system for converging Jeffery- 
Hamel flow is illustrated in Figure 4-lb. The particular geometry 
used in this study was a converging channel with a half-angle a of 
5 degrees. A polar coordinate system, centered at the intersection of 
the wall and symmetry plane, was used as the orthogonal coordinate 
system for the analysis. In Figure 4-lb, R^ is the radial coordinate 
measured from the center of the polar system, and 9 is the angular 
coordinate measured in the clockwise direction from the symmetry plane. 
Since the flow is in the -R^ direction, the computational streamwise 
coordinate is 



where R. is the radial coordinate at the initial station. The corn- 
el 

putational cross-stream coordinate is simply 
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n = f/a 

The metric scale coefficients are then 


h x = h 3 = 1 
h 2 “ R S a 

This flow, unfortunately, violates one of the basic assumptions 
made in formulating a set of equations that can be solved by spatial 
marching. The streamwise pressure gradient cannot be represented by 
an inviscid pressure gradient plus a one-dimensional viscous correction 
if self-similarity is to be maintained, for this reason the exact so- 
lution was used to compute the imposed pressure gradient P^U,n). 

The derivation of the equation for P from the exact solution for u 
is presented in Appendix H. Under these conditions, the computed vis- 
cous pressure correction p^ should be essentially zero throughout the 
duct. 

The Reynolds number used in Reference 60 was defined as 


where u^ is the centerline velocity at a given R^. For this case, 
a value of Reg = -5000 was used. (The value is negative because u CL 
is in the negative direction.) An initial channel half-wioth of 
O.Ubl cm (0.002 ft) was used as the reference length i_ r . The value 
of R<- was thus 0.6995 cm (0.02295 ft). Standard atmosphere values 
of 288 K (519’ R) and 1.2246 kg/m 3 (0.07645 lb m /f t 3 ) were used for the 
reference temperature and density, respectively. The centerline veloc- 
ity at the initial station, which for Reg = -5000 was 10.45 .../sec 
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(34.29 ft/sec), was used for the reference velocity. The reference 
Reynolds number was thus 435.8. 

The exact solution was used for the initial streamwise velocity 
profile. (See Appendix H for the equation). The initial dimensionless 
values of cross-flow velocity, density, and temperature were zero, one, 
and one, respectively. No-slip and no-bleed boundary conditions were 
used for u and v at the wall. An adiabatic wall boundary condition 
was used for the temperature. Fifty-one grid points were used between 
the symmetry plane and the wall. The points were lightly packed near 
the wall with = 10 in equation (3.1). One hundred marching 
steps were taken with a step size at of 0.05. The final area was 
thus 0.56 times the initial area. The computational grid is shown in 
Figure 4-17 with every other streamwise station plotted. The inset in 
the figure shows the grid near the wall magnified by a factor of five. 
A circle is centered at corresponding locations in the full grid and 
the inset. This calculation required 28.9 seconds of CPU time on an 
IBM 370/3033. 

4.2.2 Results 


The computed streamwise velocity profiles are presented in Figure 
4-13. Profiles are plotted every other marching station. The flow 
accelerates through the duct and the boundary layer thins somewnat. 
Self-similarity is maintained in the calculation, as shown in Figure 
4-19, where the computed solution after 100 marching steps is compared 
with the exact solution. The two solutions are the same to within the 
resolution of the plot. 
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The cross-flow velocity profiles are presented in Figure 4-20. 
Except for a slight starting transient, where it reaches a dimension- 
less value of about 0.0005, the cross-flow velocity is essentially 
zero everywhere in the duct. 

The viscous pressure correction also remains essentially zero. 

It reaches a maximum at the last station of about 0.0018 (nondimen- 
sionally), compared to an overall level of about 757. The viscous cor- 
rection to the pressure gradient at this station is -0.0017, compared 
to the exact values of -0.49 at the centerline and 0.73 at the wall. 




*v* 
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4.3 TURBULENT DEVELOPING PIPE FLOW 


One of the most basic test cases for turbulent flow is developing 
flow in a circular pipe. Although no exact solution exists, the prob- 
lem has been studied experimentally by several investigators (Refs. 47, 
61-64). The data of Barbin (Ref. 62) were chosen for comparison with 
the present analysis. 

4.3.1 Grid and Initial Conditions 


As in the laminar developing pipe flow case, a cylindrical coor- 
dinate system was used, as illustrated in Figure 4-1. Again, the com- 
putational U,n) and physical (z,r) coordinates coincide, so that 
t = z and n = r. Also, as in the laminar case, the imposed pres- 
sure P(i,n) is equal to a constant everywhere in the duct. The 
streamwise pressure gradient is completely determined by the viscous 
correction p^U). 

The pipe radius in Reference 62 was 10.16 cm (4.0 in). This was 
used as the reference length L f in the analysis. The reference 
velocity U r was 29.9 m/sec (98.1 ft/sec), the average velocity in the 
pipe. The reference temperature T p that was used was the standard 
atmosphere value of 288 K (519* R). The Reynolds number in Reference 
62, based on pipe diameter and average velocity, was 388,000. In order 

3 

to match this, the reference density was set equal to 1.143 kg/m 
(0.07138 lb m /ft^). The reference Reynolds number was thus 194,000. 

The measured velocity profile at the first experimental station, 
1.5 diameters downstream of the pipe entrance, was used as the initial 
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streamwise velocity profile for the analysis. Weighted average quadra- 
tic interpolation was used to obtain values between the tabulated ex- 
perimental data points. In addition, the method outlined in Section 
2.2.3 was used to obtain realistic velocity values between the wall 
and the first experimental point away from the wall. The initial 
cross-flow velocity v was set equal to zero. The initial dimension- 
less density and temperature were both set equal to one. This case 
was run using both the Cebeci-Smith and McDonald-Camarata turbulence 
models in the inner region. 

No-slip and no-bleed boundary conditions were used for u and v 
at the wall. An adiabatic wall boundary condition was used for the 
temperature. Fifty-one grid points were used between the centerline 
and the wall, with the points packed near the wall using equation (3.1) 
with = 100. A constant marching step size ac of 0.1 was used. 
The grid near the pipe entrance is plotted in Figure 4-21a. The grid 
was the same through the entire pipe. The boxed region near the wall 
in Figure 4-21a is shown magnified 20 times in Figure 4-21b. A circle 
is centered at corresponding locations in the two figures. 

In the experiment of Reference 62, the last data station was 
8.84 m (29.0 ft), or 87 radii, from the pipe entrance. The data indi- 
cate that the flow, although close, was not quite fully developed at 
this point. The calculations, therefore, were extended to 120 radii 
downstream of the entrance. Thus, 1170 marching steps were needed 
between 5=3 and 5 = 120. The value of A5 could probably be 
greatly increased, especially downstream of the near-entrance region, 
to cut down the number of steps, and therefore computer time, required 
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for this case. The case as run took 4.1 mi notes of CPU time on an IBM 
370/3033. 

4.3.2 Results 


The computed streamwise velocity profiles are shown in Figure 4-22. 
Profiles are plotted every 20 marching steps, with the abscissa dis- 
placed by 1/4 of the interval width for each successive profile. These 
results were computed using the Cebeci-Smith turbulence model in the 
inner region. The differences between these results and those obtained 
using the McOonald-Camarata turbulence model were too small to be dis- 
tinguished in this type of plot. The boundary layer growth along the 
wall can be clearly seen in Figure 4-22. The edge of the boundary 
layer reaches the pipe centerline about halfway through the pipe. From 
this point through to the end of the pipe the velocity profile changes 
shape very gradually. The dimensionless centerline velocity, for 
example, only increases from 1.210 at t « 70 to 1.21b at t « 120. 

The computed cross-flow velocity profiles are presented in Figure 
4-23. As in the laminar developing pipe flow case, the analysis 
quickly generates a realistic profile even though the initial 
v-velocity was set equal to zero. The point of maximum v-velocity 
moves away from the wall as the boundary layer thickness increases. 

Once the flow becomes completely viscous, the v-velocity quickly de- 
creases and remains essentially zero. Although no experimental values 
for the cross -flow velocity are given in Reference 62, these profiles 
are qualitatively the same as those obtained experimentally and ana- 
lytically by Richman and Azad (Ref. 47). 



The u-velocities computed using both the Cebeci-Smith and 
McDonald-Camarata turbulence models are compared with the data of Ref- 
erence 62 in Figures 4-24 and 4-25, respectively. Streamwise velocity 
is plotted as a function of axial distance at five radial locations. 

The agreement between theory and experiment is generally very good. 

The two turbulence models give essentially the same results, except 
very near the wall where the McDonald-Camarata model is in slightly 
better agreement with the data. These results are also compared with 
the data of Reference 62, in the form of velocity profiles at five 
axial locations, in Figures 4-26 and 4-27. 

It should be noted that in the present analysis the centerline 
velocity continually increases and approaches its fully-developed value 
asymptotically from below. Several authors (Refs. 64-67) have reported 
that the centerline velocity in developing turbulent pipe and channel 
flow typically overshoots its fully-developed value slightly, and then 
approaches it from above. This overshoot phenomenon is not apparent 
in the data of Reference 62, perhaps because measurements were not 
made far enough downstream. The exact cause of this behavior is not 
clear, but it apparently cannot be predicted with an algebraic turbu- 
lence model, at least without some strictly empirical problem- dependent 
modifications (Refs. 6b, 68). 

The computed and experimental pressure coefficients are plotted as 
a function of axial distance in Figure 4-28. Recall that since the im- 
posed pressure in this case is a constant, this pressure coefficient is 
determined by the viscous correction p'. The agreement between theory 
and experiment is fail good, with the experimental results falling 
about halfway between the computed results from the two turbulence 
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models. This again indicates that the method used to compete the vis- 
cous pressure correction is valid, and that the slight disagreement in 
Figure 4-28 is an effect of the turbulence model. 

In Figure 4-29, the computed ski» friction coefficient is plotted 
as a function of axial distance for both turbulence models. Tne shear 
stress at the wall was found by second-order one-sided numerical dif- 
ferentiation of the computed streamwise velocity. Also shown in the 
figure is the expected fully-developed value computed from Prandtl's 
formula for wall friction (Ref. 28). This formula is 


where 



A 


8t. 


F-D 



4c 


f 


(4.1) 


A curve fit to equation (4.1), accurate to ±3 percent, is given by 
White (Ref. 28) as 

a , 1.02 (log Re 0 )" 2,5 (4.2) 

Equation (4.2) is the one actually used to get the fully-developed skin 
friction value presented in Figure 4-2y. 

With both turbulence models, the skin friction decreases rapidly 
near the pipe entrance, where toe boundary layer thickness is in- 
creasing rapidly. Unlike the results in laminar developing flow, 
however, it dips below the final fully-developed value and then gradu- 
ally increases. Although no direct skin friction measurements were 
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made in the experiment of Reference 62, and the data are not detailed 

enough to compute reliable \alues, these results are qualitatively 

similar to those obtained in other experiments (Refs. 66-67). Of the 

two turbulence models, the results obtained using the Cebeci-Smith 

model are in better agreement with equation (4.2). 

Finally, the predicted Reynolds stress profiles are compared with 

data at four axial stations in Figure 4-30. The Reynolds stress, in 

2 

the form u ' v ' / (u ) , is plotted as a function of normalized dis- 
t F-U 

tance from the wall. Here u is the friction velocity for fully- 

t F-D 

developed flow, given by 



For the predicted values of Reynolds stress, equation (4.2) was used to 

compute t y . The experimental values were taken directly from Ref- 
"F-D 

erence 62, where the value of t u is not given. The experimental 

r-D 

and theoretical values presented in Figure 4-30 may therefore be norma- 
lized by slightly different values. It's interesting to note that the 
Reynolds stress curves for the two turbulence models would essentially 
coincide if each were normalized using its own predicted value of u 

t F-U 

The agreement between theory and experiment is qualitative, at 
best, at z/R « 9 and 33, but is fairly good at z/R ■ 57 and 81. In 
the experiment, transition from laminar to turbulent flow was artifi- 
cially promoted using a strip of sand grain roughness near the pipe 
entrance. This could be the cause of the disagreement at the first 
two stations in Figure 4-30. Of course, the disagreement could also 
be due simply to inadequacies in the turbulence models. 
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•4.4 S—SHAPEP DUCT FLOW 


The next test case studied was turbulent flow in a two-dimensional 
S-shaped duct. This is a fairly complex flow, with the boundary layer 
on each wall seeing both favorable and adverse pressure gradients. In 
addition, the curvature terms in the equations become important since 
the metric scale coefficients vary in both coordinate directions. The 
particular geometry used for this case v as the 30*-45* S-duct studied 
experimentally by Butz (Ref. 69), 

4.4.1 Geometry arid Computational Coordinates 

The geometric configuration is presented in Figure 4-31. Two 
circular arc bends, both with a centerline radius of 0.508 m (20.0 in), 
were used to make up the S-duct. The first bend covered an arc of 30*, 
and the second an arc of 45*. The duct width was 0.1016 m (4.0 in). 

Tne computational coordinate system for this configuration was 
computed using the method of References 2 arid 3, as described in Sec- 
tion 2.1. In order to satisfy the uniform inflow and outflow boundary 
conditions that are used with this method, straight sections four 
inches, or one duct width, long were added upstream and downstream of 
the S-duct. The resulting coordinate system is shown in Figure 4-32. 
Distances have been nondimensionalized by the outer wall radius. (As 
a check, a coordinate system was later constructed using straight sec- 
tions that were two duct widths long. The metric parameters for the 
two cases were essentially identical). The coordinates were generated 
with 51 cross-stream points and 100 streamwise points, evenly spaced 
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(in the computational coordinates { and n) In both directions. 

During the viscous marching solution, the metric parameters at each 
point In the viscous grid Mere found by Interpolation over this evenly 
spaced grid. 

The metric coefficients generated by the method of References 2 
and 3 are, unfortunately, not smooth, in Figure 4-33 the metric scale 
coefficient h (recall h^ - h ? with this method) is plotted along 
both walls as a function of axial distance. It was found that the com- 
puted streamwise velocities were not affected by the oscillations shown 
in Figure 4-33a. These oscillations were, howevet, reflected In the 
computed cross-flow velocities, although the calculation did not go un- 
stable. To eliminate the sharp oscillations in the metric coeffi- 
cients, they were smoothed by splitting the h - h(c) curve at each 
constant n-line into several sections, separated by "knots," and using 
a cubic spline smoothing within each section. The knot locations were 
chosen to yield tl° best overall fit, in the least squares sense, to 
the original h - h(t) curve at n - 0. The same knot locations 
were then used for all th» n-lines. (New knot locations could not be 
used for each n-line because this led to sharp changes in h in the 
n-direction at a given t). The smoothed metric coefficients, using 
10 knots, are plotted in Figure 4-33b. The smoothing removed the sharp 
oscillations, although a couple of stretched out "bumps" remain. These 
caused no difficulties in the viscous calculation, however. 
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4.4.2 Imposed Pressure Field 

Since the computational mesh consists of the streamlines and 
potential lines from a planar potential flow solution, the imposed 
pressure was computed directly from the metric scale coefficients, as 
described in Section 2.3.2. The imposed, or inviscid, pressure coef- 
fic.ent on each wall is compared with the experimental data in Figure 
4-34. Here c p is defined as 

P - P re f 

c p ' r — .T" 

7 “A 

where p p is the density, U r is the average velocity at the duct 
entrance, and p ref is the static pressure at the first measurement 
location on the upper wall. 

The gradients given by the inviscid pressure distribution match 
the experimental pressure gradients very well. In particular, the 
elliptic effects of the geometry on the pressure field are correctly 
reproduced. The pressure along the top wall starts to decrease ahead 
of the inflection station, at 30 degrees, and increases ahead ot the 
exit, at 75 degrees. The pressure along the bottom wall, of course, 
does just the opposite. This imposed pressure field, therefore, cor- 
rectly introduces the elliptic effects of the geometry into the viscous 
marching solution. 

One additional point should be mentioned here. In the experiment 
of Reference 69, the S-duct was constructed from two separate circular 
arc sections. Flow surveys were made by inserting a 3.81 cm (1.5 in) 
long survey section, containing a translating probe, at either end of 
a circular arc section. Since straight sections were added upstream 
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and downstream of the S-duct when generating the computational mesh 
(and therefore the pressure field), the effect of this survey section 
on the flow is properly accounted for when comparing with the data at 
the entrance and exit of the S-duct. Putting the survey section at 
the Inflection station, however, changes the geometry slightly. This 
change turns out to be important, as will be shown later, because the 
Inflection station is where the streamwise pressure gradients are the 
largest. A second computational mesh and Imposed pressure field were 
therefore generated, with the survey section between the two circular 
arc sections, for use when comparing with the data at the inflection 
station. 

4.4.3 Initial Conditions 


The reference conditions used in the analysis were chosen to 
match the flow conditions in the experiment of Reference 69. The 
average inlet velocity of 19.2 m/sec (63 ft/sec) was used as the 
reference velocity U r . The reference length L f was 0.5588 m 
(1.83333 ft), the outer wall radius at the inlet. To match the exper- 
imental conditions, a reference Reynolds number Re f of 7.216x10^ was 
then required. The reference temperature was assumed to be 288 K 

*3 

(519* R), and the reference density was computed to be 1.2042 kg/nr 
(0.075173 lb m /ft 3 ). 

The measured streamwise velocity profile at the entrance station 
was used to start the marching procedure. This profile is shown in 
figure 4-35. Interpolation was used to get values between the data 
points. The initial cross-flow velocity was set equal to zero, ana 


the initial dimensionless density and temperature Mere both set equal 
to one. The marching calculation was started 3.33 cm (1.312b in) up- 
stream of the start of the first bend. This is the actual measurement 
location for the entrance station, accounting for the presence of the 
survey section and the 1.429 cm (0.5625 in) upstream extension of the 
probe itself (Ref. 70). 

No-slip and no-bleed boundary conditions were used for the veloc- 
ities at the walls. An adiabatic wall boundary condition was used for 
the temperature. Fifty-one grid points were used between the two 
walls. Because the boundary layers were thin, the points were tightly 
packed near both walls, with ■ 200 in equation (3.1). The stream- 
wise step size at was 0.05. The actual mesh used in the viscous cal- 
culation, for the case without the survey section in the middle, is 
shown in Figure 4-36. For clarity, only every other streamwise station 
is shown. The inset in the figure shows the grid near the wall magni- 
fied by a factor of 2b. A circle is centered at corresponding loca- 
tions in the full grid and the inset. For this case 163 marching steps 
were needed to reach the end of the duct, requiring 4a. b seconds of CPU 
time on an IBM 370/3033. For the case with the survey section, 170 
marchiny stops were needed and 50.6 seconds of CPU time were required. 

4.4.4 R esults 

This case was run using both the Cebeci-Smith and McOonald- 
Canarata turbulence models. The results were virtually identical, so 
only the results obtained for the Cebeci-bmith turbulence model will 
be shown in this section. 



lot) 

The computed streamwise velocity protiles are shown in figure 
4-3/a. The boundary layer growth through the duct can be clearly seen* 
although it's still relatively thin at the exit station. The slope of 
the velocity profile in the essentially inviscid core region is consis- 
tent with a potential vortex flow, with the velocity inversely propor- 
tional to the distance from the center of curvature of the duct sec- 
tion. Ihis can be seen more clearly in Figures 4-3/b and 4-3/c, where 
the predicted profiles midway through the first and second bends are 
plotted to a larger scale. The velocity near the bottom wall (i.e., 
the inner wall in the first bend and the outer wall in the second beno) 
increases as the flow enters the first bend, decreases around the in- 
flection station, and increases again leaving the second bend. Near 
the top wall, of course, the opposite occurs. 

The computed streamwise velocity profiles at the inflection and 
.nit stations are compared with the experimental profiles in Figure 
4-38. The computed values at the inflection station are for a config- 
uration with the survey section between the two circular sections, as 
described in Section 4,4,2. The actual locations of the inflection 
ana exit stations, taking into account the survey section and the up- 
stream extension of the probe, are 0.47o cm (0.18/S in) downstream of 
the ends ot the first and second bends. In general the agreement be- 
tween theory and experiment is very good. 

It is interesting to note how important it is to include the ef- 
fects ot the survey section and probe, especial 1; when comparing with 
the data at the inflection station. The velocity profile changes 
shape rapidly in this region because of the steep favorable anu ad- 
verse pressure gradients on the top and bottom walls, respectively. 
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In Figure 4-39, the theoretical results exactly at the inflection 
plane, ignoring the effects of the survey section, are compared with 
the data. The agreement is not nearly as good as that shown in Figure 
4-38a. 

In Figure 4-40 the pressure coefficient for the vi ous flow along 
each wall is compared with the experimental data. The pressure coeffi- 
cient is determined by the imposed, or inviscid, pressure plus the 
viscous pressure correction. By comparing with Figure 4-34, it can be 
seen that the viscous correction in the first bend yields slightly 
worse agreement with the data along the top wall and about the same 
level of agreement along the bottom wall. In the second bend the 
agreement is improved along the top wall and is slightly worse along 
the bottom wall. It should be noted, however, that because of the 
scatter in the data, it is difficult to determine exactly how well the 
theory matches the experiment. In defining the pressure coefficient, 
the pressure at the first measurement location on tne top wall was 
used as the reference pressure. If this particular measured pressure 
were too high or low, all of the experimental values would be too low 
or high relative to the predicted values. 
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Figure 4-H Imposod pressure cooffidonts compared «4th <Ma for S-thopod dud ftoa. 



S-shapod dud flow. 
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4.5 SUBSONIC DIFFUSER FLOW 

As a final test case, compressible turbulent flow in an axisym- 
metric subsonic diffuser typical of those used in supersonic propulsion 
systems was studied. The diffuser chosen is called a dM/dz diffuser 
because the area distribution is such that the one-dimensional Mach 
number distribution through the diffuser is linear. This type of dif- 
fuser has been tested as a component of a supersonic inlet system 
(Kef. 71), as part of an investigation of various diffuser types 
(Ref. 6), and as part of a vortex generator study (Kef. 72). Because 
of their yreater detail, the data of Shaw (Ref. 72) were chosen for 
comparison with the present analysis. Two cases were run, both with 
the Cebeci-Smith turbulence model, with average inlet Mach numbers c 
0.29 and 0.75. 

4.5.1 Geometry and Computational Coordinates 

The geometric configuration is presented in Figure 4-41. The cowl 
radius was constant at 15.24 cm (6.0 in), while the centerboay radius 
varied from 11.43 cm (4.5 in) at the diffuser inlet station, where 
z = 0.0, to 5.387 cm (2.121 in) at the exit station, where z = 22.9 
cm (9.0 in). Upstream of the diffuser itself was a slightly diverging 
throat section, with a centerbody wall angle of 3 degrees. The diffu- 
ser area ratio was 2.0. In this type of diffuser the area change, and 
thus the rate of diffusion, is fairly small near the entrance, but 
increases rapidly in the downstream half of the diffuser. Also shown 
in Fiyure 4-41 are the locations of th« five total pressure rakes. 
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designated A to £, that were used in the experiment of Keference 72. 
These rakes were vertical, and were located at i « 0.0, b.90, 11.4, 

17. «, and 22.9 cm (0.0, 2.7b, 4.b, 7.0, and 9.0 in). 

The computational coordinate system was again computed using the 
method described in Section 2.1. A cubic fairing was added along the 
centerbody in the region of the diffuser exit to eliminate the sharp 
change in wall curvature at that point, which would have caused numer- 
ical problems in generating the coordinate system. Straight sections 
were again added upstream and downstream of the duct to satisfy the 
uniform inflow and outflow boundary conditions assumed by the coordi- 
nate generation method. A longer straight section was used at the down- 
stream end than at the upstream end because of tne larger annular width 
and the larger centerbody curvature at that point. The resulting coor- 
dinate system is shown in Figure 4-42. Distances have been nondimen- 
sionalized by the cowl radius. Near the diffuser exit plane, the ac- 
tual centerbody contour is drawn as a dashed line, showing the extent 
of the cubic fairing used in this region. The coordinates were gener- 
ated using 51 cross-stream points and 100 streamwise points, evenly 
spaced (in the computational coordinates £ and n) in both direc- 
tions. During the viscous marching solution, the metric parameters at 
each point in the viscous grid were found by interpolation over this 
evenly spaced grid. 

The metric scale coefficients along each wall are shown in Figure 
4-43. Along the straight outer cowl the distribution is smooth. Along 
the curved centerbody, however, there are oscillations. It was again 
found, as in the S-duct case of Section 4.4, that the computed stream- 
wise velocities were not affected by these oscillations. Therefore, 
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tor the results presented in this section* no smoothing of the metric 
coefficients was performed. 

4.5.2 imposed Pressure field 

The imposed, or inviscid, pressure field was computed using the 
method described in Section 2.3.2. Since this is an axisymmetric case, 
this involved solving the ax i symmetric potential equation, equation 
(2.38), to get an incompressible inviscid velocity. The lieblein- 
Stockman correction was used to include the effects of compressibility 
on t tie imposed pressure field. The resulting pressure coefficients 
along each wall are compared with the experimental data, for both the 
* 0.2y and = 0.75 cases, in figures 4-44 and 4-45, respectively. 
Here c p is defined as ( P-P ret )/Up where is the nominal inlet 
dynamic pressure and p refr tor each wall is the measured static pres- 
sure closest to the diffuser inlet station. For the cowl this was 
exactly at the irlet station and for the centerbody it was 2.54 cm 
(l.u in) upstream. 

In general the agreement between the inviscid pressure coeffi- 
cients and the data is fairly good. I he bumpiness in the pressure 
distribution along the centerbody is a result of the bumpiness in the 
metric scale coefficient, lhe viscous pressure correction computed 
during the marching solution will tend to shift the theoretical values 
down as axial distance increases. Note the flattening of the experi- 
mental pressures along the centerbody at the downstream end tor both 
Mach numbers. ibis indicates the presence of boundary layer separation 
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on the centerbody near z/L - 1.1. No separation is indicated along 
the cowl . 

4.i>.3 initial Conditions 


The reference conditions used in the analysis were chosen to match 
the flow conditions in the experiment of Reference 72. The cowl radius 
was used for the reference length L r . The average inlet values of 
velocity, temperature, and density were used for the reference values. 
These reference conditions are summarized in Table 4-1 for the two 
cases. 

The marching calculation was started at the diffuser inlet, the 
location of rake A in the experiment. The initial boundary layer pro- 
file could not be determined in detail from the data, however. There- 
fore the initial streamwise velocity profile was constructed by speci- 
fying the free-stream velocity plus boundary layer thickness parameters 
on both the cowl and centerbody. The velocity in the bounaary layers 
was then found using the procedure for turbulent flow described in 
Section 2.2.3. The boundary layer parameters were adjusted by judi- 
cious trial-and-error to give a good fit to the experimental profile. 
The values used are presented in Table 4-2. 

The initial cross-flow velocity was set equal to zero. The ini- 
tial dimensionless density and temperature were both set equal to one 
in the tree-stream. The values of density and temperature in the 
boundary layer were found using the procedure of section 2.2.3. 

No-slip and no-bleed boundary conditions were used for the veloc- 
ities at the walls. An adiabatic wall boundary condition was used tor 
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the temperature. Fifty-one grid points were used In the cross-stream 
direction, packed at both walls with 0^ - SO in equation (3.1). 

The streamwise step size at was 0.0S. The actual mesh used in the 
viscous calculation is shown in Figure 4-46. The inset in the figure 
shows the grid near the wall magnified by a factor of 25. A circle is 
centered at corresponding locations in the full grid and the inset. 

For both Mach numbers the analysis predicted flow separation before 
the end of the diffuser was reached. For the ■ 0.29 case, 87 

marching steps were taken, requiring 27.0 seconds of CPU time on an 
IBM 370/3033. For the M^ . 0.75 case, 78 marching steps and 24.2 
seconds of CPU time were used. 

4.5.4 Results 


The computed streamwise velocity profiles are shown in Figures 
4-47a and b. For both inlet Mach numbers, the boundary layers on the 
cowl and centerbody grow rapidly in the downstream half of the diffuser 
until the flow finally separates along the centerbody. The predicted 
separation location, taken as the point where negative streamwise ve- 
locities first occur, was z/L p • 1.26 for the * U.29 case 
and z /L r - 1.10 for the M^ - 0.7b case. Both of these locations 
are in good agreement with those indicated by the experimental pressure 
distributions shown in Figures 4-44a and 4-45a. The inarching procedure 
was able to continue past separation, but it quickly became unstable, 
giving physically unrealistic results. The results shown in Figure 
4—47 • therefore, include profiles only up to the separation point. It 
should be noted that by modifying the streamwise convective terms in 
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the governing equations a marching procedure can be made stable in 
separated flow regions (Ref. 73). However, the resulting equations 
are approximate and can be used only if the separation region is small 
(i.e., the flow quickly reattaches) and does not appreciably affect 
the rest of the flow field. This is clearly not the case for the flows 
considered here. 

The computed streamwise velocity profiles for the - 0.29 
case are compared with the experimental results in Figure 4-48. The 
velocity ratio is plotted as a function of fractional distance across 
the duct. Comparisons can only ce made at the first four rakes since 
the flow separates before the last rake is reached in both the analysis 
and the experiment. The agreement between the predicted and experimen- 
tal.results is generally very good. 

The analytical results for the Mj - 0.7b case are compared 
with the data in Figure 4-49. Since the analysis predicted separation 
slightly upstream of the rake U location for this case, the predicted 
results can only be shown for the first three rakes. The agreement 
with the data at these rakes is again very good. The experimental ve- 
locity profile at rake 0, which is also shown in Figure 4-49, does 
not indicate separated flow. However, the experimental static pressure 
distribution along the centerbody (see Figure 4-45a) indicates separa- 
tion close to the rake D location. It's possible that the flow ac- 
tually was separated at rake u, but that the t irst tube in the rake 
was too tar from the wall to detect it. 

The pressure coefficient for the viscous flow along each surface 
is compared with the experimental data in Figure 4-50 for the 
Mj * 0.29 case ami in Figure 4-bl for the Mj = 0.7b case, by 
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comparing with Figure 4-44, It can be seen that the viscous blockage 
correction shifts the pressure downward only slightly for the ■ 
0.29 case. There is a larger shift In the ■ 0.75 case (compare 
Figures 4-51 and 4-45) because of the thicker boundary layer along the 
centerbody. The analytical and experimental pressure gradients along 
the centerbody match closely, but the absolute pressure levels them- 
selves are slightly off. This could be due to not having exactly the 
right reference pressure when computing the pressure coefficients. 
Along the cowl, the predicted and experimental pressures agree well 
near the beginning of the diffuser, but disagree further downstream. 
This is probably due to the upstream influence of the separated flow 
region in the experiment. The large separation region reduces the 
effective flow area, causing the flow along the cowl to accelerate. 
This effect cannot be predicted, of course, by a marching procedure. 

Finally, the theoretical and experimental Reynolds stresses at 
rakes B, C, and D for the - 0.29 case are shown in Figure 

4-52. (The Reynolds stresses were not. measured fcr the « 0.75 
case, or at rake A for the - 0.29 case). From the figure, it 

can be seen that the Cebeci-Smith turbulence model yields at least 
qualitatively correct Reynolds stresses at rakes B and C. The 
agreement is much poorer at rake D, which is very near the point of 
boundary layer separation. 
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Table 4-1. Reference Conditions for Diffuser Flow Test Cases 


1 

i Mi « 0.29 

i 

i i 

i Mi - 0.75 i 

i i 

i i 

i L r i 0.152 m (0.5 ft) 

I 1 

i i 

i 0.152 m (0.5 ft) i 

i U r i 99.40 m/sec (326.1 ft/sec) 

i 247.0 m/sec (810.3 ft/sec) i 

i p r i 2.0055 kg/m 3 (0.12520 lb m /ft 3 ) 

i | 

i 1.6087 kg/m 3 (0.10043 lb m /ft 3 ) i 

i 1 

i T r i 292.4 K (526.3* R) 

i 269.9 K (485.8* R) i 

i i 

i Re r i 1.68xl0 6 

i i 

) 1 

i 3.56xl0 6 i 

i i 


Table 4-2. Initial Profile Parameters 
for Diffuser Flow Test Cases 



Ml * 0.29 

i Mi = 0.75 

i 

(u/U r ) e 

1.05 

i 

i 1.03 

i 

( 6/ L r )C0WL 

0.03 

i 0.03 

i 

( 6*/ L r )coWL 

0.00333 

i 0.00429 

I 

(®/ L r)c0WL 

0.00267 

1 

i 0.00321 
| 

(«/L r )c-B0DY 

0.0433 

i 0.0433 
| 

(«*/L r )c-B0DY 

0.00394 

i 0.00867 
| 

( e / L r )c— BODY 

0.00328 

» 0.00578 
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(a) Cantarbody. 

Flgura 4-44. - Imposad pra**ura coefficients eomparad with data far 
diffuser flow, Mj ■ CL29. 
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(a) Csfltarbody. 
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Figure 4-45. -Conclude 


126 


ORIGINAL PAGE IS 
OF POOR QUALITY 



*i/j ‘sunt* 


taMMna.tn, 

F^urt 4-46. - CoMpuWteml «rM u«4 for viscous cotoutatton. Offusor ftar. 


127 


ORIGINAL p 
0F POOR qi 


it 





Fractional dhtonca across dud Fractional Astanct across dud 


128 


ORIGINAL PAC2 13 
OF POOR QUALITY 



0 .2 .4 .6 .1 LO LO LO LO 

ValocKy ratio, uAJ r 

Figuro 4-41 - Computed and nparlmantal velocity profiles for dlffusar flow, Mj • a 29. 



0 .2 .4 .4 .1 LO LO LO 

Velocity ratio, uAJ r 

Figure 4*40. • Computed and axptrlmantel velocity profiles ter dlffusar floor, Mj • & 7i 


LO 


> *>>t> > 


Pressure cttfffctont, Cp Pressure cwffteltnt, Cp 



Pressure coefflcterrt, cu 



131 



SECTION 5 


CONCLUSIONS 

A new solution procedure nas been developed and used to compute 
compressible viscous subsonic flow in both planar and axisymmetric 
ducts. A set of equc-' «cns that can be solved by forward marching was 
derived by neglecting second derivatives in the streamwise direction 
and by uncoup 1 ing the streamwise and cross-flow pressure gradients. 
The streamwise pressure gradient was written as the sum of a known two- 
dimensional imposed pressure gradient and an unknown one-dimensional 
viscous correction computed as part of the marching procedure. The 
governing equations were solved simultaneously using an implicit 
finite-difference method, based on the results of the test cases pre- 
sented in Section 4, a number of specific conclusions can be made about 
the analysis. 

1. The analytical results agree extremely well with experimental 
data and/or exact solutions for the laminar flow cases studied. This 
indicates that the basic solution procedure is valid since the turbu- 
lence model, of course, has no effect in these cases. 

2. The analytical results for turbulent flow also agree well with 
data, although not quite as well as for laminar flow. Both turbulence 
models that were used give equally good results for the cases studied. 
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3. The method used to compute the viscous pressure gradient cor- 
rection is valid. This was demonstrated by the laminar and turbulent 
pipe flow cases, where the streamwise gradient is given solely by the 
viscous correction, and by the Jeffery-Hamel flow case, where the vis- 
cous correction is essentially zero. 

4. When the imposed pressure field is a good approximation to the 
actual pressure field, the computed results will be accurate. This 
was demonstrated by the Jeffery-Hamel flow case. 

5. The analysis can be used for flows with both favorable and ad- 
verse pressure gradients. It can also be used to accurately predict 
the location of flow separ'tion, as shown by the diffuser flow test 
cases. 

b. The predicted streamwise velocities are not especially sensi- 
tive to any oscillations present in the metric scale coefficients and 
imposed pressure field. These oscillations are reflected in the pre- 
dicted cross-flow velocities, but not to the point of numerical in- 
stabi 1 ity. 

/. The marching procedure can be started with the initial cross- 
flow velocity set equal to zero when the actual profile is unknown, as 
is often the case in engineering practice. Ihe analysis will then 
generate a realistic profile in the first few marching steps. 

8. The analysis is fast enough for use in parametric studies and 
design work. A typical turbulent flow case with a blxlUO mesh will run 
in about 30 seconds on an IBM 370/ 3033 . 
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APPENDIX A 


DERIVATION OF GOVERNING EQUATIONS 


The general equations for compressible* viscous fluid flow, in 
orthogonal curvilinear coordinates, can be found in one form or 
another in several references (e.g., References 74-75). Here they are 
written as: 


CONTINUITY 
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In these equations, h^ .hg, and h^ are the metric scale coeffi- 
cients for the orthogonal curvilinear coordinate system; £, n, and e 
are the three coordinate directions; u, v, and w are the velocities 
in the t, n, and $ directions, respectively; p is the static 
density; p is the static pressure; and h is the static enthalpy. 

The subscripts c, n» e, and t denote partial differentiation. 

The shear stresses are given by 



and is the molecular viscosity. The second coefficient of vis- 
cosity, x , 's assumed equal to -2w, / 3. 
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The heat fluxes are given by 


«i - - ^ k i T < 
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^ k L T n 




(A.7) 


where T is static temperature and k ( is the molecular thermal 
conductivity. 

Finally, the viscous dissipation, , is given by 


> L - M L jz (e^ + + e 33 ) ♦ Ue ?3 ) 2 + Ue 31 )‘ 


+ ^ e 12^] “ 3 \ (e ll + e 22 + e 33^ 
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where 
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Equations (A.1)-(A.5) are, in theory, valid tor turbulent as well 
as laminar flow. In practice, however, solving them directly for tur- 
bulent flow would be a hopeless task. The extremely fine resolution 
that would be required in both space and time is simply beyond the 
capacity ot current and foreseeable computers. 

The way around this problem is to replace the instantaneous quan- 
tities in the equations with the sum of their average and fluctuating 
values. This can be aone using either conventional time averaged 
values or mass averaged values. The conventional time average of a 
quantity t is defined by 
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f dt 


(A. 9) 


where b is some time period that is long compared to the character- 
istic time scale of the turbulent fluctuations. The mass averaged 
value of a quantity f is defined by 


f = 




(A. 10) 


Here, following the procedure of Cebeci and Smith (Ref. 43), mass 
averaged values are used for the velocities and enthalpy, and conven- 
tional time averages are used for the remaining variables. Therefore, 


u = u + u' 
v * v + v' 
w = w + w' 


h = h + h* y 

T = T + T' 

p = p + p" 


(A. 11) 


J 

Note that an overbar and double prime have been used for conventional 
time averaging and a tilde and single prime for mass weighted averag- 
ing. With mass weighted averaging, as shown in Reference 43, 


p7 = pf 

pf‘ = 0 l 

V = 7 - f 


(A.12) 
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The expressions for u, v, w, h, p, and p are substituted into 
equations (A.1)-(A.5), and the resulting equations are then time 
averaged. At this point, the equations are restricted to steady, 
two-dimensional, mean flow. The results are as follows: 
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(A. 16) 


ENERGY 



Because of the inherently three-dimensional nature of turbulence, 
there are correlation terms involving w‘ appearing in these two- 
dimensional mean flow equations. In fact, at this point, the 
s-momentum equation has not been completely eliminated. 

Following the experimental and order-of-magnitude arguments of 
Reference 43, several assumptions are now made. 

1. Fluctuations in pressure, viscosity, thermal conductivity, and 
specific heat are neglected. 

2. Triple correlations are neglected. This gives 



pu v 


rrrr _ ~ 


etc. 
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3. Ti lues can be replaced by overbars on u, h, and T, but not 
on v. Thus, 


u • u 
h ■ K 
T ■ T 


p 

Assuming T - T also means that the time averaged equation of state 
for a perfect gas, given by 


P * R pT 

can be written as 

p - R p T (A. 18) 

The energy equation can be modified by using the definition of 
enthalpy, 

h « e ♦ -E- (A. 19) 

p 

— ae ae 

where e is specific internal energy, and then replacing — ana — 

at an 

3 I a y 

with c v — and c v — . Note that this does not imply constant c^, 
only a perfect gas. 

lhe iPlp r and v‘p terms can be eliminated using the momentum 
t n 

equations. When this is done, several turbulence terms appear that are 
analagous to the terms m the viscous dissipation, f, (Ret. 7b). To- 
gether, these terms are called the apparent dissipation, 
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The continuity equation can be used to eliminate the p ceriva- 
tives in the energy equation and to write the momentum and energy 
equations in nonconservative form. When all of this is done, equa- 
tions {A. 13)— {A- 17) become: 

CONTINUITY 
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It is now assumed that the Reynolds stresses (p u‘v* , etc.) can be 
modeled in exactly the same form as the viscous stresses (t^. etc.) 

in equation (A. 6), by using a turbulent viscosity p-j.. Thus, 


- p u ' 2 - 2 »*t(^ u 5 + 7“! W T V ' V 
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etc. The turbulent viscosity is defined by 
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where e is an eddy viscosity which comes from some appropriate tur- 
bulence model. Note that this formulation assumes the eddy viscosity 
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is an isotropic scalar. At this point, the e-momentum equation drops 
out for two-dimensional flow. It is also assumed that the turbulent 


heat fluxes (p FF and p v'h* ) in the energy equation can be modeled 
in the same form as the molecular heat fluxes (q^ and q^ ) in equa- 
tion (A. 7) by using a turbulent thermal conductivity, ky. Thus, 


where 
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When these modeling assumptions are incorporated into equations 
(A. 20 )— (A.22) and (A. 24), and the overbars and tildes dropped for sim- 
plicity, the results are as follows: 
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APPENDIX B 

ORDER-OF -MAGNITUDE ANALYSIS 

In order to use a forward marching solution procedure, second de- 
rivatives in the t direction, which appear in the shear stress 

| 

terms, must be eliminated. An order-of-magnitude analysis is used to j 

show that this is justified. First, it is assumed that the flow is f 

% 

primarily in the t direction. Thus v, the velocity in the n | 

3 

i 

-t 

direction, is assumed to be of order 6. (Here 6 is taken to be a j 

small number representing the order-of-magnitude difference between 
the streamwise and cross-flow velocities. It is not necessarily a J 

boundary layer thickness, since, in this internal flow analysis, vis- . 

cous terms are included throughout the flowfield and the boundary lay- 
ers can gr w to fill the duct.) It is expected that gradients normal 
to the walls will be greater than gradients in the streamwise direc- | 

tion. Therefore, derivatives in the n direction are assumed to be I 

of order 1/6. In order to make the largest viscous term the same order 

as the convective terms, the effective viscosity is set equal 

to order 6 . (Equivalently, the Reynolds number is assumed to be 

of order 1/6 .) Similarly, k E is also set equal to order 6 . 

All other terms, including the metric coefficients and their deriva- 
tives, are set equal to order 1. 
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With these assumptions, the magnitude of the effective shear 
stresses can be estimated as shown below. The order of magnitude is 
written below each term. 
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Then, in the streamwise momentum equation, 
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In the cross-flow momentum equation. 
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And in the energy equation, 
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All terms of order 1 and order 6, at least, are to be retained. 
Therefore, in the streamwise momentum equation, the term ^2 h 3 T E ) 
can be eliminated. In the cross-flow momentum equation, the 0(« 3 ) 
part of the term ^ h 2 h 3 T E ) can *^ so be e, * m i nated * is 

enough to eliminate all second derivatives with respect to t, and 
thus allow the use of a marching solution procedure. Many of the 
other terms are high order and could also be eliminated, but it is not 
necessary since they can be handled within the framework of a numeri- 
cal marching procedure. 
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APPENDIX C 


POTENTIAL FLOW COMPRESSIBILITY CORRECTION 


For compressible flow, the incompressible velocity given by equa- 
tion (2.40) is modified using the Lieblein-Stockman compressibility 
correction (Ref. 37). This correction was developed specifically for 
internal flows. The compressible velocity is given by 


V 


C 



(C.l) 


where V. is the average incompressible velocity at a o 
’AVE 

station, p, is the incompressible (or total) density, and p r is 

U AV£ 

the average compressible density at a station. One-dimensional 

isentropic relations are used to get p-/ p r . Thus, 

1 L AVE 




l/u-l) 


(C.2) 


where a n is the total speed of sound. The value of V r /a n is 
u C AVE u 


found from 


a* h * n (l,n)/C2(lrl)3 c ave 
r ■ (V-) 


a o L 1 *^ 


fif)‘ 


1/(t-1) 


(C.3) 


where 
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O 

** poor 


p°cj is 

Walit? 



Y / , v(Y*l)/(T-lf 

1/2 

.* (H J 

Po 


Equation (C.3) is solved for V r /a n by Newton iteration. 

l AVE u 


APPENDIX 0 


COUPLED DIFFERENCE EQUATIONS 

In this appendix, the governing equations are written in finite- 
difference form. These equations are the result of the differencing 
ana linearization procedures presented in Sections 3.1 and 3.2. The 
equations correspond to the differential equations ( 2 -13)— ( 2.16) , with 
the shear stress and heat flux terms written out in full. Since the 
equations in this appendix are valid at any j-location in the grid, 
the "j" subscript has been omitted. In acu lion, the metric scale 
coefficients are unaerstood to be at the i + w,j) grid location, so 
these subscripts are also omitted. The a notation from Section 3.2 
is used, so that ap * p. + j - p., etc. The variable in the 

equations is given by 

u i+w ' "i * "‘("i - "i-l> 

The subscripts t ano n denote partial differentiation. 
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By expanding the a and 6^ operators, and after much collect- 
ing of terms, these equations can be written in the following forms. 
The dependent variables in equations (l).5)-(D.8) are at the unknown 
i+1 station, and the M i+1" subscript has been omitted. 
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As described in Section 3.2, the effective viscous dissipation i 
evaluated at station i and thus appears in the source term in 
the energy equation. Note also that the source term S x in the 
streamwise momentum equation contains the viscous pressure correction 
computed during the preliminary marching step described in Sections 
2.3.1 and 3.5. Although the “j" subscript has been omitted from the 
coefficients and source terms, they will be different at each 
j-location in the grid. 

Equations (D.5)-(D.8) are written for j = 2 to j = J-l. The 
values of the dependent variables at j = 1 and j = J are either 
known functions of 5 or expressible in terms of the dependent 
variables at other j-locations (see Section 3.3). These equations 
thus represent a system of 4(J-2) coupled linear algebraic equations 
in 4 ( J— 2 ) unknowns. They can also be written as a matrix equation 
with a block tri-diagonal coefficient matrix. 
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APPENDIX E 

DENSITY BOUNDARY CONDITION AT A WALL 

The density at a solid surface is found by evaluating the 
cross-flow momentum equation at the surface. The difference equation 
at the wall is the same as at interior points, and is presented in 
Appendix D as equation (D.3). At the wall, however, the 6y and 
6y operators in equation (D.3) are given by equations (3.16) and 
(3.18) . When the resulting algebraic equation is solved for p^ , 
it has the form 



Vj±i 


+ C 


W, p j±2 


* Vj + 


C W 6 u j±l 


+ c, u 

*6 


W, u j±2 


+ C 


w 7 v j 


+ c, , v + L, v, i0 + C,, T, 


W 8 J±1 % T j±2 


10 


w n ^ 


c w 


12 


J±2 


+ S, 


(£. 1 ) 


where j - 1 or J. 

I 

The C w s and can be expressed in terms of the 

I 

C y s and $ y presented in Appendix D. Their exact form 
depends on the type of temperature boundary condition used. For a 
specified wall temperature. 
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C u « C v i 3C V 

W 1 '2 Y i 


Cy ■ t 4Cy 

"2 T 1 


Cu * + Cw 
w 3 T 1 


C - - C Y * 3C y 
w 4 5 t 4 


Cy — ± 4Cy 
"5 4 


Cy = + Cy 

w 6 4 


.r ? 3C ♦ J_ 4 1 Wl , M „ , F , ,> 

\ S Re r^UJ 1 ‘ 2 ' 


± 4C 


• 2 4 1 {M\ w ^ 


1 r , + 1 4 1 /d¥\ 2 „ it 

S W “ Ui TZ7 


Cy - ~ Cy + 3 C y 
w 10 T 11 t 9 


Cy - t 4 C y 
W 11 T 9 


Cy * ^ Cy 
"12 9 


S W “ S Y 
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where 


r. r .14 1 / dY y w AC 

S ■ % 3 7 Z w Ui UV) 2 

In these equations the top sign is used for the j ■ 1 boundary and 
the bottom sign is used for the j ■ J boundary. 

For a specified temperature gradient normal to the wall, C w ^, 


C y -C u , and S u change. They become 
"10 "12 


\ - C Y 2 - ""‘Vi* At * wrt \ (w) w + w 


"i+1 





AC ± 3Cy 


c w = 0 

w i0 


4 r 

’“n ’ * 7 Cy n 


c w = \ C Y 

"12 J 11 


( E .2b) 


S W ' S Y * 7 (lr) u (sV)C 


y "11 

i + l 
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APPENDIX F 


COEFFICIENT MATRICES AND SOURCE TERMS 


When the boundary conditions presented in Section 3.3 are applied 
to equations (D.5MD.8), the resulting equations can be written as 


B 2 x 2 ♦ c 2 X 3 . S 2 


(F.la) 


Vj-i * 5 o V Vi ’ o < j < o - 2) (F.ib) 


A J-1 X J-2 + B J-1 X J-1 ‘ S 0-’ 


(F.lc) 


where the A's, B's, and C's are the known 4x4 coefficient matrices, 
the S's are the known 4-element source term vectors, and the X's are 
the 4-elemer: solution vectors given by 


p 

u 

v 

T 


(F.2) 


For 3 <. j <. 0-2, the coefficient matrices and source term vectors are 
given by 


1 


1 


10 


10 


10 


10 


(F • 3a) 
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The elements in the above matrices are defined in Appendix D. 

The form of the coefficient matrices and source terrr vectors at 
j = 2 and j = J-l depend on the type of boundary conditions being used. 
In the following discussion, let (B^ represent the m'th row of B^. 

This notation will also be used for , A j_i» and c » 



18S 


ORIGINAL Pt\QZ r~ 
OF POOR QUALITY 


For a symmetry line at j « 1, the boundary conditions given in 

equations (3*19) are applied to equations (D.5)-(D.8). Each row of 
the matrix then has the form 



(F.4a) 


For convenience, (B 2 ) m has been written as the transpose of a column 
vector. In this equation, C represents C^, Cy, and C £ for 
m = 1, 2, 3, and 4, respectively. This convention will be followed 
throughout this appendix. Each row of C 2 has the form 



C - 4 C 

eq 3 3 eq x 

c - — c 

eq 6 3 U eq 4 



eq 


12 


1 r 

1 So 


(F .4b) 


The source term vector S 2 is unchanged from the form presented in 


equation (F.3d). 
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At a wall with a specified temperature, the boundary conditions 
given by equations (3.20), (3.21), (3.22a), and (3.23) are applied to 
equations (D.5)-(D.8). For a solid surface at j • 1 with a specified 
temperature, then, each row of B£ has the form 






4 8 c 

V eq i 



(F . 5a) 


and each row of has the form 


V. 


c c 

eq 3 ^ eq l 


C ♦ C 

eq 6 ^ eq l 


C + ^ C 
e q 9 eqj 


"12 

C + T r ~~~ 

eq 12 V eq l 


(F.5b) 
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Similarly, each element of the source term vector has the form 




1 * ( e "4 + S C *«i) Vl ' (s 


V'(So*V c,qi ) T " 1 


For a specified wall temperature at j « J, each row of Aj_^ has the 


r + - r 

S x eq 3 


(Aj_i) m - 


r + o p 

eq 4 ^ eq 3 


C ^7 + X Ceq3 


"12 

C + — c 
e< >10 C Wj e «3 
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and each row of has the form 





♦ 




C ♦ C 
eq 5 eq 3 

C W 

c ♦J!§C 
e <8 eq 3 




(F.6b) 


Each element of Sj ^ has the form 


‘ViL 


eq 


w r 


r + -_Z. r 
eq 6 ^ eq 3 


- C 




eq g 


'7 


- C 


^ eq^y \ et lv> c u e q *j i w 




TO 


T2 "W 


1 




(F.6c) 


At a wall with a specified temperature gradient normal to the 
wall, the boundary conditions are given by equations (3.20), (3.21), 
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(3.22b), and (3.23). For a specified temperature gradient at the 
j * 1 wall then, each row of § 2 has the form 

C, 


< S 2>m * 


U 


c _ + 


^ e0 i 
Cw 5 

C * q 5 + ^ S 

r 

Ce, a * S S 


” e «n + C ^i + * C ^ioj 


(F. 7a) 


and each row of C 2 has the form 


* 


W 3 

r + - .. c 

eq 3 C Wi eq 


1 


W 6 

c + ^ c 

eq 6 ^ eq l 


w g 

C e, 9 * ^ C e 9l 


w 

c ♦ 12 c 1 c 

eq 12 ^ 6q l 1 eq 10 


(F.7b) 
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Each element of has the form 


(S,L - S 


w 


'2'm * 5 eq “ eq 1 


K + ^ Ceq iJ V 


-I c 


eq 7 


+ ^s)vr * 


10 


(ty) w lY 

\ > W i+1 


(F.7c) 


For a specified temperature gradient at the j « J wall, each row 
of has the form 




W 3 

e <l ^ et >3 


W 6 

C e, 4 * C eq 3 


C._ * -r^-C 


eq 7 eq 3 

C W 

+ 12 C 1 C 

eq 10 \ eq 3 1 eq 12 


(F . 8a) 
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and each row of Bj_^ has the form 


S * ^ S 


(B J-l J m 


C + yr~— C 
eq 5 eq 3 


C n + C 
eq 8 \ eq 3 


C + C + d C 
eq ll C W X eq 3 ^ eq 12 


Each -lenient of Sj^ has the form 


(S J-l ) m " S eq “ C C eq 3 “l C eq 6 + C C *%)\ n ' ( C eq g 


s >3/ Vi 


When a wall function boundary condition is used for the stream- 
wise velocity, and the temperature gradient normal to the wall is 



192 


ORIGINAL PAGE 13 
OF POOR QUALITY 


specified, the boundary conditions are given by equations (3.21), 
(3.22b), (3.23), and (3.24). Under these conditions at the j - 1 
wall, each row of EL still has the form given in equation (F.7a). 

^ «w 

The form of each row of C£ becomes 



W 3 

C ♦ ft — — C 

eq 3 eq x 

C W. + C W fi 

C + C + — ^ - C 

eq 4 eq 6 ^ eq l 

C W 

"g 

C ♦ 7C- 2 - C 

eq, eqj 


C .Jllic -Ac 

eq 12 ^ eq l 7 e<) 10 


<F.9») 


In addition, the form of each element of becomes 



S e, * ^ \ 



W 4 

+ -p — ■ c 
\ eq li 








aY 


+ 


(F .9b) 
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At the J - 0 wall, each row of still has the form given in 

equation (F.8b). The form of each row of becomes 




W 3 

C ♦ c 
e «l \ eq 3 

Gu + 

W 4 W 6 

C + C ♦ C 

eo 4 eq 6 ^ eq 3 


S + ^ C *<3 


'eq 


10 


"12 1 

eq 3 3 eq 12 


(F • 10a) 


and the form of each element of Sj ^ becomes 


(S J-l ) m 85 S eq “ C eq 3 ' (^ C eq 6 + C eq 3 J F u w ’ ^eq g 


'W 


7 


^ C e q 3 y Vw i+i ^ Ceq 12 V/ W i+1 


aY 


(F .10b) 


The parameter F appearing in equations (F -9b) and (F.lOb) is 


defined by 


F u - (2 aY) 
U W 


(aro/ 


i*i T j y o+i 


(F.ll) 
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where the top sign is used for the j - 1 wall and the bottom sign is 
used for the j - J wall. 

When a wall function boundary condition is used for the stream- 
wise velocity, and the wall temperature is to be specified, the 
boundary conditions are given by equations (3.21), (3.23) (3.24) , and 

(3.25). Under these conditions at the j ■ 1 wall, each row of Bg 
has the form 
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and each row of Cg has the form 



In addition, the form of each element of becomes 







+ 


^7 

^ eq i 




+ 




Vt. 


(F . 12b) 


(F .12c) 
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For these boundary conditions at the j ■ 0 wall, each row of Aj j 
has the form 



and each row of has the form 
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The form of each element of Sj_^ becomes 



The parameter F T is defined by 

‘w 



(F.13c) 


(F.14) 
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APPENDIX G 

DIFFERENCE EQUATIONS FOR VISCOUS PRESSURE CORRECTION 

In this appendix, the difference equations used to compute the 
viscous pressure correction p* during the preliminary marching step 
are derived. First the streamwise momentum equation, presented in 
Section 2.2.1 as equation (2.14), is uncoupled from the continuity, 
cross-flow momentum, and energy equations and written in finite- 
difference form. The uncoupling is done by evaluating p, u, v, and 
T at the known station i. The "j" subscript has been omitted from I 

all terms in the equation, since it is valid at any j-location. In I 

1 

addition, the "i+w" subscript has been omitted from the metric scale f 

coefficients. The notation from Section 3.2 is used, so that 

* * 

au » U 'j+}~ u i* The superscript * is used on u to distinguish 

it from the streamwise velocity to be computed during the main march- 
ing step. Note, however, that the variables at the known station i 
are those from the previous completely-coupled main marching step. 

The subscripts t and n denote partial differentiation. The 
resulting uncoupled streamwise momentum difference equation is as 

I 

follows: j 

! 
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h 3 h l 
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_1_ u + 1_ dY 
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Continued on next page 
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* 

u-, - u 


h 2. h 3 


Tl l* *2 ' *3 ) +W AU) f 


x A- 

T3 


2 1 


u^, - Uj 


h 3 h 2 


T W 1 — “ * I KT ( 2 K? " E* 


'1 \ "3 "2 



After collecting terms and applying the boundary conditions, 
equation (G.l) can be written as 


h 2 u Z + c 2 u 3 * d 2 p C ‘ S X. 


a j u j-i * b j u j * c jVi + d j p i * s Xj < 3 i j i J - 2) (G - 2) 


®J-l u J-2 + b J-l u 0-l + d J-l p i " S X 


J-l 


where, for 2 <_ j J-l, 
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Now, for 2 <_ j < J-l» lot 
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1 2 1 dY 7 
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1 dY 
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1 C 4 1_ dY 
Re r h 1^2 3 h ° dn Vi 
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3h 2 lh 1 h 3 K v i 


h 2 h 3 


+21 i . 2 1 I 0 K C \ x 

3 h x u i AC 3 hj_ l 2 h 2 ~ h 3 ) {i " w V u i Li 


h 3 r u / h 3 h ? 
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Continued on next page 
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2 1 dY . u + 2 1 ( 7 

Vi *iiq l 2 



At 


For 3 < j <_ 0-2, the coefficients in equation (G.2) are simply 
given by 

a . « a ' 

J J 




C j “ C j 


S X = s x 

x j j 


However, the forms of the coefficients b^, c^, a j_j» an( ^ ^j-l* anc * 
the source terms s ¥ and s Y depend on the type of boundary 

2 J - 1 

condition being used. For a symmetry line at j = 1, using equation 
(3.19) gives 


b 2 ■ b 2 * 3 a 2 
C 2 * C 2 " 1 a 2 
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For a solid wall at j - 1 with u^ - u 


W* 


b 2 «b* 


C 2“ C 2 


S X 2 “ S X 2 ~ a 2 u W 


For a wall function boundary condition at j * 1, using equation 
(3.24) gives 


b 2- b 2 


C 2 = a 2 + c 2 


S X, - s x 2 + *2 \ 


where is defined by equation (F.ll). For Uj ■ u w at j * 0, 

U W 


a J-l = a J-l 


b J-l ‘ b J-l 


X 0-l Vl ^ 


And finally, for a wall function boundary condition at j = J, 


a J-l - a 0-l + C J-1 


b J-l * b J-l 


S V * Sy + c' . F, 
X J-1 X J-1 0-1 °W 
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The other equation used tc compute the viscous pressure correc- 
tion Is the total mass flow rate equation, derived in Section 3.5.2. 
It is repeated here. 



(G.3) 


•k 

The unknowns in this equation are u, + , and p' . By moving the 

1 *i+w 

known terms to the right side, equation (G.3) can be rewritten as 




Equation (6.4) is to be integrated numerically from j ■ 1 to 
j = J using Simpson's rule. Since Simpson's rule requires an odd 
number of grid points, if J is even the integration from 0-1 to 
J is done by trapezoidal rule. When the integration is performed, 
and the boundary conditions are applied, the resulting equation can be 
written as 

e 2 u 2 + e 3 u 3 + •** * e 0-2 u J-2 + e J-l u J-l + fp e “ S M 
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where 


' • ^ ■’ * W " 


( 2 ») 


Now, for 1 <_ j <, J» let 


ei « 


j-KT7( P M + >i (S3?),*1 


s,n * (2*)“ ay 




Then for even j's from 4 to J-3, the coefficients in equation (G.5) 
are given by 


and for odd j's from 5 to J-3, 


e j “ 4e j 


- 2e i 


The exact forms of e^ and depend on the type of 
boundary condition used at j = 1. For a symmetry line at j 


1. 


4 

e-, * 4eA + ->r e- 

.. _ L. 


2 “ "~2 1 1 
3 


e, * 2e^ - y e x 


For a solid surface at j * 1 with u^ = u^. 


e 2 " 4e 2 


e 3 * 2e 3 


upttMftwi «*>**' IV *•"" 'I' 1 * 


207 


PAGE © 

F POOR QUALITY 


And for a wall function boundary condition at j - 1, 


e 2 ■ 4e 2 


e, - e{ ♦ 2e$ 


The forms of ej_ 2 and ej_., depend on the boundary 

condition used at j * J, and on whether J is odd or even. For 
J odd and Uj = u w , 

e J-2 ‘ 2e j-2 


e J-l * 4e j-l 


For J odd with a wall function boundary condition at j ■ J, 

e 0-2 = 2e G-2 * e J 


e J-l " 4e J-l 


For J even and Uj * u^ , 


e J-2 * 4e J-2 


e J-l = f e j-l 


For J even with a wall function boundary condition at j ■ J, 


a — 4p 1 + — p 1 

e J-2 " 4e 0-2 2 e J 


e j-i * I e j-i 


The form of the source term s M depends on the types of 
boundary conditions used at j « 1 and at j * J, and on whether 
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J is odd or even. First, for J odd with symmetry conditions at 
j - 1, 


s m “ s m " e J F J 


v 

where Fj - u^ for a specified velocity boundary condition at j - J, 

I 

and F. - -F for a wall function boundary condition at j » J. For 
J Uj 

J odd with a solid surface at j ■ 1, 

S M “ S M ‘ e i F i “ 

I 

where F^ « u^ for a specified velocity boundary condition at 

l 

j = 1, and F, ■ - F for a wall function boundary condition at j » 1. 
i u l 

For J even with a symmetry line at j - 1, 

5 m * s m ' I e 3 F j 

And finally, for J even with a solid surface at j ■ 1, 

s m - S M " e i F i - I e o F i 
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APPENDIX H 


JEFFERY-HAMEL FLOW STREAMWISE PRESSURE GRADIENT 


One of the classical examples of viscous duct flow Is laminar 
Incompressible flow In a two-dimensional wedge-shaped channel, known 
as Jeffery-Hamel flow. When self-similarity is assumed, an exact 
solution to the Navier-Stokes equations exists for this flow. The 
streamline' are radial, intersecting at a line source (for diverging 
flow) or sink (for converging flow). The solution for the velocity, 
derived in terms of Jacobian elliptic functions, is given by Mill saps 
and Pohlhausen (Ref. 60). From Reference 60, the velocity profile for 
converging Jeffery-Hamel flow is given by 

F(») - 2 [m 2 (k 2 - 2) - 1] «■ 6 m 2 (l - k 2 ) dn“ 2 (mp,k) (H.l) 


where 



(H.2) 


m 


2 


1 + Re 0 /2 

n- 

1 - 2k* 


(H.3) 


and k is the largest number satisfying the equation 


Re 0 (l - 2 k J ) 

2k ? ^3k ? - 3 ♦ ^ (ic 2 _ 2)j 


2 

sn (mo.k) 


(H.4) 
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Here sn(ma,k) and dn(m*,k) are Jacobian elliptic functions, which 
are evaluated by the Arithmetic-Geometric Mean method of Reference 
77. The Reynolds number ReQ is given by 



u CL R $ 


The polar coordinates R $ and 9 are defined in Figure 4-16, 
Uq L is the centerline velocity at a given R s , and a is the 
half-angle of the duct. For the case studied here (Rep ■ -5000, 
o * 5 degrees) the value of k is 0.99988175. 

The radial momentum equation for purely radial flow, in the 
coordinate system of Figure 4-16, is 


u 



1 3D . / 3^U , 1 3U 

v « 


♦ 


1 3^U U \ 


( H .5 ) 


Substituting equation (H.2) for u gives 



(H.6) 


Using derivatives of Jacobian elliptic functions as given by Refer- 
ence 77, the second derivative on the right-hand side of equa- 
tion (H.6) is given by 


« 12m 4 (l - k 2 )k dn~ 2 (m*,k) [- sn 2 (m*,k) dn 2 (m*,k) 
69 L 


+ ci 2 (m*,k) dn 2 (mf,k) ♦ 3k sn 2 (m*,k) cn 2 (m*,k)J 


211 


Here cn(m*,k) Is also a Jacobian elliptic function. 
Since the computational streamwise coordinate t 


C 



R 


S 


the exact streamwise pressure gradient is 


ifi. 

H 




is given by 


(H.7) 


